Context Navigation

source: codes/icosagcm/trunk/src/write_field.f90 @ 169

Last change on this file since 169 was 161, checked in by dubos, 11 years ago
Minor additions
File size: 62.7 KB

Line
1	module write_field
2	USE genmod
3	implicit none
4
5	INTEGER,SAVE :: ncprec
6
7	TYPE ncvar
8	INTEGER :: size
9	INTEGER,POINTER :: nc_id(:)
10	INTEGER :: displ
11	END TYPE ncvar
12
13	integer, parameter :: MaxWriteField = 1000
14	integer, dimension(MaxWriteField),save :: FieldId
15	TYPE(ncvar), dimension(MaxWriteField),save :: FieldVarId
16	integer, dimension(MaxWriteField),save :: FieldIndex
17	character(len=255), dimension(MaxWriteField) :: FieldName
18
19	integer,save :: NbField = 0
20
21	contains
22
23	SUBROUTINE init_writeField
24	USE ioipsl
25	use netcdf_mod
26	IMPLICIT NONE
27	CHARACTER(LEN=255) :: netcdf_prec
28
29	netcdf_prec='float'
30	CALL getin("netcdf_prec",netcdf_prec)
31
32	SELECT CASE(TRIM(netcdf_prec))
33	CASE('float')
34	ncprec=NF90_FLOAT
35	CASE('double')
36	ncprec=NF90_DOUBLE
37	CASE default
38	PRINT*,'Bad selector for variable netcdf_prec : <', TRIM(netcdf_prec),"> options are <float>, <double>"
39	STOP
40	END SELECT
41	END SUBROUTINE init_writeField
42
43	function GetFieldIndex(name)
44	implicit none
45	integer :: GetFieldindex
46	character(len=*) :: name
47
48	character(len=255) :: TrueName
49	integer :: i
50
51
52	TrueName=TRIM(ADJUSTL(name))
53
54	GetFieldIndex=-1
55	do i=1,NbField
56	if (TrueName==FieldName(i)) then
57	GetFieldIndex=i
58	exit
59	endif
60	enddo
61	end function GetFieldIndex
62
63	SUBROUTINE Writefield(name_in,field,nind,once)
64	USE domain_mod
65	USE field_mod
66	USE transfert_mpi_mod
67	USE dimensions
68	USE mpipara
69	IMPLICIT NONE
70	CHARACTER(LEN=*),INTENT(IN) :: name_in
71	TYPE(t_field),POINTER :: field(:)
72	INTEGER,OPTIONAL,INTENT(IN) :: nind
73	LOGICAL,OPTIONAL,INTENT(IN) :: once
74	LOGICAL :: once_
75
76	TYPE(t_field),POINTER :: field_glo(:)
77
78	!$OMP MASTER
79	IF(PRESENT(once)) THEN
80	once_=once
81	ELSE
82	once_=.FALSE.
83	END IF
84
85	IF (field(1)%ndim==2) THEN
86	CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type)
87	ELSE IF (field(1)%ndim==3) THEN
88	CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3)
89	ELSE IF (field(1)%ndim==4) THEN
90	CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3,field(1)%dim4)
91	ENDIF
92
93	CALL gather_field(field,field_glo)
94
95	IF (mpi_rank==0) THEN
96	IF (PRESENT(nind)) THEN
97	CALL writefield_gen(name_in,field_glo,domain_glo,nind,once=once_)
98	ELSE
99	CALL writefield_gen(name_in,field_glo,domain_glo,1,ndomain_glo,once=once_)
100	ENDIF
101	ENDIF
102
103	CALL deallocate_field_glo(field_glo)
104	!$OMP END MASTER
105
106	END SUBROUTINE writefield
107
108	! SUBROUTINE Writefield(name_in,field,nind)
109	! USE netcdf
110	! USE domain_mod
111	! use field_mod
112	! USE dimensions
113	! USE geometry
114	! IMPLICIT NONE
115	! CHARACTER(LEN=*),INTENT(IN) :: name_in
116	! TYPE(t_field),POINTER :: field(:)
117	! INTEGER,OPTIONAL,INTENT(IN) :: nind
118	! REAL(r8),ALLOCATABLE :: field_val2d(:)
119	! REAL(r8),ALLOCATABLE :: field_val3d(:,:)
120	! REAL(r8),ALLOCATABLE :: field_val4d(:,:,:)
121	! TYPE(t_domain),POINTER :: d
122	! INTEGER :: Index
123	! INTEGER :: ind,i,j,k,n,ncell,q
124	! INTEGER :: iie,jje,iin,jjn
125	! INTEGER :: status
126	! CHARACTER(len=255) :: name
127	! CHARACTER(len=255) :: str_ind
128	! INTEGER :: ind_b,ind_e
129	! INTEGER :: halo_size
130	! LOGICAL :: single
131	!
132	!
133	! name=TRIM(ADJUSTL(name_in))
134
135	! IF (PRESENT(nind)) THEN
136	! name=TRIM(name)//"_"//TRIM(int2str(nind))
137	! PRINT *,"NAME",nind,int2str(nind),name
138	! ind_b=nind
139	! ind_e=nind
140	! halo_size=1
141	! single=.TRUE.
142	! ELSE
143	! ind_b=1
144	! ind_e=ndomain
145	! halo_size=0
146	! single=.FALSE.
147	! ENDIF
148
149	! Index=GetFieldIndex(name)
150	! if (Index==-1) then
151	! call create_header(name,field,nind)
152	! Index=GetFieldIndex(name)
153	! else
154	! FieldIndex(Index)=FieldIndex(Index)+1.
155	! endif
156	!
157	! IF (Field(ind_b)%field_type==field_T) THEN
158	! ncell=1
159	! DO ind=ind_b,ind_e
160	! d=>domain(ind)
161	! IF (Field(ind)%field_type/=field_T) THEN
162	! PRINT *,"Writefield, grille non geree"
163	! RETURN
164	! ENDIF
165
166	! n=0
167	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
168	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
169	! IF (d%own(i,j) .OR. single) n=n+1
170	! ENDDO
171	! ENDDO
172
173	! IF (field(ind)%ndim==2) THEN
174	! ALLOCATE(Field_val2d(n))
175	! n=0
176	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
177	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
178	! k=d%iim*(j-1)+i
179	! IF (d%own(i,j) .OR. single) THEN
180	! n=n+1
181	! Field_val2d(n)=field(ind)%rval2d(k)
182	! ENDIF
183	! ENDDO
184	! ENDDO
185	! status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d, &
186	! start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /))
187	! DEALLOCATE(field_val2d)
188	! ELSE IF (field(ind)%ndim==3) THEN
189	! ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
190	! n=0
191	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
192	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
193	! k=d%iim*(j-1)+i
194	! IF (d%own(i,j) .OR. single) THEN
195	! n=n+1
196	! Field_val3d(n,:)=field(ind)%rval3d(k,:)
197	! ENDIF
198	! ENDDO
199	! ENDDO
200	! status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /), &
201	! count=(/n,size(field(1)%rval3d,2),1 /))
202	! DEALLOCATE(field_val3d)
203	! ELSE IF (field(1)%ndim==4) THEN
204
205	! DO q=1,FieldVarId(index)%size
206	!
207	! ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
208	! n=0
209	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
210	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
211	! k=d%iim*(j-1)+i
212	! IF (d%own(i,j) .OR. single) THEN
213	! n=n+1
214	! Field_val3d(n,:)=field(ind)%rval4d(k,:,q)
215	! ENDIF
216	! ENDDO
217	! ENDDO
218	! status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /), &
219	! count=(/n,size(field(1)%rval4d,2),1 /))
220	! DEALLOCATE(field_val3d)
221	! ENDDO
222	! ENDIF
223	!
224	! PRINT *,NF90_STRERROR(status)
225	! ncell=ncell+n
226
227	! ENDDO
228	!
229	! ELSE IF (Field(ind_b)%field_type==field_Z) THEN
230	! ncell=1
231	! n=0
232	! DO ind=ind_b,ind_e
233	! d=>domain(ind)
234	! CALL swap_geometry(ind)
235	! CALL swap_dimensions(ind)
236	!
237	! n=0
238	! DO j=jj_begin+1,jj_end
239	! DO i=ii_begin,ii_end-1
240	! n=n+1
241	! ENDDO
242	! ENDDO
243
244	! DO j=jj_begin,jj_end-1
245	! DO i=ii_begin+1,ii_end
246	! n=n+1
247	! ENDDO
248	! ENDDO
249
250	! IF (field(ind)%ndim==2) THEN
251	! ALLOCATE(Field_val2d(n))
252
253	! n=0
254	! DO j=jj_begin+1,jj_end
255	! DO i=ii_begin,ii_end-1
256	! n=n+1
257	! k=iim*(j-1)+i
258	! Field_val2d(n)=field(ind)%rval2d(k+z_down)
259	! ENDDO
260	! ENDDO
261
262	! DO j=jj_begin,jj_end-1
263	! DO i=ii_begin+1,ii_end
264	! n=n+1
265	! k=iim*(j-1)+i
266	! Field_val2d(n)=field(ind)%rval2d(k+z_up)
267	! ENDDO
268	! ENDDO
269
270	! status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1), &
271	! Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /))
272	! DEALLOCATE(field_val2d)
273
274	! ELSE IF (field(ind)%ndim==3) THEN
275	! ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
276	! n=0
277	! DO j=jj_begin+1,jj_end
278	! DO i=ii_begin,ii_end-1
279	! n=n+1
280	! k=iim*(j-1)+i
281	! Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:)
282	! ENDDO
283	! ENDDO
284
285	! DO j=jj_begin,jj_end-1
286	! DO i=ii_begin+1,ii_end
287	! n=n+1
288	! k=iim*(j-1)+i
289	! Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:)
290	! ENDDO
291	! ENDDO
292	! status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /), &
293	! count=(/n,size(field(1)%rval3d,2),1 /))
294	! DEALLOCATE(field_val3d)
295	! ELSE IF (field(1)%ndim==4) THEN
296
297	! DO q=1,FieldVarId(index)%size
298	! ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
299	! n=0
300	! DO j=jj_begin+1,jj_end
301	! DO i=ii_begin,ii_end-1
302	! n=n+1
303	! k=iim*(j-1)+i
304	! Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q)
305	! ENDDO
306	! ENDDO
307
308	! DO j=jj_begin,jj_end-1
309	! DO i=ii_begin+1,ii_end
310	! n=n+1
311	! k=iim*(j-1)+i
312	! Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q)
313	! ENDDO
314	! ENDDO
315
316	! status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,1,FieldIndex(Index) /), &
317	! count=(/n,size(field(1)%rval4d,2),1 /))
318	! DEALLOCATE(field_val3d)
319	! ENDDO
320	! ENDIF
321	!
322	! PRINT *,NF90_STRERROR(status)
323	! ncell=ncell+n
324
325	! ENDDO
326	!
327	! ENDIF
328	! status=NF90_SYNC(FieldId(Index))
329	!
330	! END SUBROUTINE Writefield
331
332
333	SUBROUTINE Writefield_gen(name_in, field, domain_type, ind_b_in, ind_e_in,once )
334	USE netcdf_mod
335	USE domain_mod
336	USE field_mod
337	IMPLICIT NONE
338	CHARACTER(LEN=*),INTENT(IN) :: name_in
339	TYPE(t_field), POINTER :: field(:)
340	TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:)
341	INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in
342	INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in
343	REAL(r8),ALLOCATABLE :: field_val2d(:)
344	REAL(r8),ALLOCATABLE :: field_val3d(:,:)
345	REAL(r8),ALLOCATABLE :: field_val4d(:,:,:)
346	LOGICAL, INTENT(IN) :: once
347	TYPE(t_domain),POINTER :: d
348	INTEGER :: Index
349	INTEGER :: ind,i,j,k,n,ncell,q
350	INTEGER :: iie,jje,iin,jjn
351	INTEGER :: status
352	CHARACTER(len=255) :: name
353	CHARACTER(len=255) :: str_ind
354	INTEGER :: ind_b,ind_e
355	INTEGER :: halo_size
356	LOGICAL :: single
357
358	name=TRIM(ADJUSTL(name_in))
359
360	IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN
361	name=TRIM(name)//"_"//TRIM(int2str(ind_b))
362	ind_b=ind_b_in
363	ind_e=ind_b_in
364	halo_size=1
365	single=.TRUE.
366	ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
367	name=TRIM(name)//"_"//TRIM(int2str(ind_e))
368	ind_b=ind_e_in
369	ind_e=ind_e_in
370	halo_size=1
371	single=.TRUE.
372	ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
373	ind_b=ind_b_in
374	ind_e=ind_e_in
375	halo_size=0
376	single=.FALSE.
377	ELSE
378	halo_size=0
379	ind_b=1
380	ind_e=ndomain
381	single=.FALSE.
382	ENDIF
383
384	Index=GetFieldIndex(name)
385	if (Index==-1) then
386	call create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once)
387	Index=GetFieldIndex(name)
388	else
389	FieldIndex(Index)=FieldIndex(Index)+1.
390	endif
391
392	IF (Field(ind_b)%field_type==field_T) THEN
393
394	ncell=0
395	DO ind=ind_b,ind_e
396	d=>domain_type(ind)
397	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
398	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
399	IF (d%assign_domain(i,j)==ind .OR. single) ncell=ncell+1
400	ENDDO
401	ENDDO
402	ENDDO
403
404	IF (field(1)%ndim==2) THEN
405	ALLOCATE(Field_val2d(ncell))
406	n=0
407	DO ind=ind_b,ind_e
408	d=>domain_type(ind)
409	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
410	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
411	k=d%iim*(j-1)+i
412	IF (d%assign_domain(i,j)==ind .OR. single) THEN
413	n=n+1
414	Field_val2d(n)=field(ind)%rval2d(k)
415	ENDIF
416	ENDDO
417	ENDDO
418	ENDDO
419
420	IF (once) THEN
421	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d, &
422	start=(/ 1 /),count=(/ncell /))
423	ELSE
424	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d, &
425	start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /))
426	ENDIF
427
428	DEALLOCATE(field_val2d)
429	ELSE IF (field(1)%ndim==3) THEN
430	ALLOCATE(Field_val3d(ncell,size(field(1)%rval3d,2)))
431	n=0
432	DO ind=ind_b,ind_e
433	d=>domain_type(ind)
434	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
435	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
436	k=d%iim*(j-1)+i
437	IF (d%assign_domain(i,j)==ind .OR. single) THEN
438	n=n+1
439	Field_val3d(n,:)=field(ind)%rval3d(k,:)
440	ENDIF
441	ENDDO
442	ENDDO
443	ENDDO
444
445	PRINT *, 'Writefield ', TRIM(name), MAXVAL(Field_val3d(:,1)), MINVAL(Field_val3d(:,1)) ! FIXME
446
447	IF (once) THEN
448	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1 /), &
449	count=(/ncell,size(field(1)%rval3d,2) /))
450	ELSE
451	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1,FieldIndex(Index) /), &
452	count=(/ncell,size(field(1)%rval3d,2),1 /))
453	ENDIF
454
455	DEALLOCATE(field_val3d)
456	ELSE IF (field(1)%ndim==4) THEN
457
458	DO q=1,FieldVarId(index)%size
459
460	ALLOCATE(Field_val3d(ncell,size(field(1)%rval4d,2)))
461	n=0
462	DO ind=ind_b,ind_e
463	d=>domain_type(ind)
464	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
465	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
466	k=d%iim*(j-1)+i
467	IF (d%assign_domain(i,j)==ind .OR. single) THEN
468	n=n+1
469	Field_val3d(n,:)=field(ind)%rval4d(k,:,q)
470	ENDIF
471	ENDDO
472	ENDDO
473	ENDDO
474
475	IF (once) THEN
476	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1 /), &
477	count=(/ncell,size(field(1)%rval4d,2),1 /))
478	ELSE
479	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,FieldIndex(Index) /), &
480	count=(/ncell,size(field(1)%rval4d,2),1 /))
481	ENDIF
482	DEALLOCATE(field_val3d)
483	ENDDO
484	ENDIF
485
486
487	ELSE IF (Field(ind_b)%field_type==field_Z) THEN
488
489	ncell=0
490	DO ind=ind_b,ind_e
491	d=>domain_type(ind)
492
493	DO j=d%jj_begin+1,d%jj_end
494	DO i=d%ii_begin,d%ii_end-1
495	ncell=ncell+1
496	ENDDO
497	ENDDO
498
499	DO j=d%jj_begin,d%jj_end-1
500	DO i=d%ii_begin+1,d%ii_end
501	ncell=ncell+1
502	ENDDO
503	ENDDO
504	ENDDO
505
506	IF (field(1)%ndim==2) THEN
507	ALLOCATE(Field_val2d(ncell))
508
509	n=0
510
511	DO ind=ind_b,ind_e
512	d=>domain_type(ind)
513	DO j=d%jj_begin+1,d%jj_end
514	DO i=d%ii_begin,d%ii_end-1
515	n=n+1
516	k=d%iim*(j-1)+i
517	Field_val2d(n)=field(ind)%rval2d(k+d%z_down)
518	ENDDO
519	ENDDO
520
521	DO j=d%jj_begin,d%jj_end-1
522	DO i=d%ii_begin+1,d%ii_end
523	n=n+1
524	k=d%iim*(j-1)+i
525	Field_val2d(n)=field(ind)%rval2d(k+d%z_up)
526	ENDDO
527	ENDDO
528	ENDDO
529
530	IF (once) THEN
531	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1), &
532	Field_val2d,start=(/ 1 /),count=(/ncell /))
533	ELSE
534	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1), &
535	Field_val2d,start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /))
536	ENDIF
537	DEALLOCATE(field_val2d)
538
539	ELSE IF (field(1)%ndim==3) THEN
540	ALLOCATE(Field_val3d(ncell,size(field(1)%rval3d,2)))
541	n=0
542	DO ind=ind_b,ind_e
543	d=>domain_type(ind)
544	DO j=d%jj_begin+1,d%jj_end
545	DO i=d%ii_begin,d%ii_end-1
546	n=n+1
547	k=d%iim*(j-1)+i
548	Field_val3d(n,:)=field(ind)%rval3d(k+d%z_down,:)
549	ENDDO
550	ENDDO
551
552	DO j=d%jj_begin,d%jj_end-1
553	DO i=d%ii_begin+1,d%ii_end
554	n=n+1
555	k=d%iim*(j-1)+i
556	Field_val3d(n,:)=field(ind)%rval3d(k+d%z_up,:)
557	ENDDO
558	ENDDO
559	ENDDO
560
561	IF (once) THEN
562	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1 /), &
563	count=(/ncell,size(field(1)%rval3d,2) /))
564	ELSE
565	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1,FieldIndex(Index) /), &
566	count=(/ncell,size(field(1)%rval3d,2),1 /))
567	ENDIF
568	DEALLOCATE(field_val3d)
569
570	ELSE IF (field(1)%ndim==4) THEN
571
572	DO q=1,FieldVarId(index)%size
573	ALLOCATE(Field_val3d(ncell,size(field(1)%rval4d,2)))
574	n=0
575	DO ind=ind_b,ind_e
576	d=>domain_type(ind)
577	DO j=d%jj_begin+1,d%jj_end
578	DO i=d%ii_begin,d%ii_end-1
579	n=n+1
580	k=d%iim*(j-1)+i
581	Field_val3d(n,:)=field(ind)%rval4d(k+d%z_down,:,q)
582	ENDDO
583	ENDDO
584
585	DO j=d%jj_begin,d%jj_end-1
586	DO i=d%ii_begin+1,d%ii_end
587	n=n+1
588	k=d%iim*(j-1)+i
589	Field_val3d(n,:)=field(ind)%rval4d(k+d%z_up,:,q)
590	ENDDO
591	ENDDO
592	ENDDO
593
594	IF (once) THEN
595	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,1 /), &
596	count=(/ncell,size(field(1)%rval4d,2) /))
597	ELSE
598	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,1,FieldIndex(Index) /), &
599	count=(/ncell,size(field(1)%rval4d,2),1 /))
600	ENDIF
601	DEALLOCATE(field_val3d)
602	ENDDO
603	ENDIF
604
605	ENDIF
606	status=NF90_SYNC(FieldId(Index))
607
608	END SUBROUTINE Writefield_gen
609
610
611
612	SUBROUTINE Writefield_mpi(name_in,field,nind)
613	USE netcdf_mod
614	USE domain_mod
615	use field_mod
616	USE dimensions
617	USE geometry
618	IMPLICIT NONE
619	CHARACTER(LEN=*),INTENT(IN) :: name_in
620	TYPE(t_field),POINTER :: field(:)
621	INTEGER,OPTIONAL,INTENT(IN) :: nind
622	REAL(r8),ALLOCATABLE :: field_val2d(:)
623	REAL(r8),ALLOCATABLE :: field_val3d(:,:)
624	REAL(r8),ALLOCATABLE :: field_val4d(:,:,:)
625	TYPE(t_domain),POINTER :: d
626	INTEGER :: Index
627	INTEGER :: ind,i,j,l,k,n,ncell,q
628	INTEGER :: iie,jje,iin,jjn
629	INTEGER :: status
630	CHARACTER(len=255) :: name
631	CHARACTER(len=255) :: str_ind
632	INTEGER :: ind_b,ind_e
633	INTEGER :: halo_size
634	LOGICAL :: single
635	INTEGER :: displ
636
637
638	name=TRIM(ADJUSTL(name_in))
639
640	IF (PRESENT(nind)) THEN
641	name=TRIM(name)//"_"//TRIM(int2str(nind))
642	PRINT *,"NAME",nind,int2str(nind),name
643	ind_b=nind
644	ind_e=nind
645	halo_size=1
646	single=.TRUE.
647	ELSE
648	ind_b=1
649	ind_e=ndomain
650	halo_size=0
651	single=.FALSE.
652	ENDIF
653
654	Index=GetFieldIndex(name)
655	if (Index==-1) then
656	call create_header_mpi(name,field,nind)
657	Index=GetFieldIndex(name)
658	else
659	FieldIndex(Index)=FieldIndex(Index)+1.
660	endif
661
662	IF (Field(ind_b)%field_type==field_T) THEN
663	ncell=1
664	DO ind=ind_b,ind_e
665	d=>domain(ind)
666	IF (Field(ind)%field_type/=field_T) THEN
667	PRINT *,"Writefield, grille non geree"
668	RETURN
669	ENDIF
670
671	n=0
672	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
673	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
674	IF (d%own(i,j) .OR. single) n=n+1
675	ENDDO
676	ENDDO
677
678	displ=FieldVarId(index)%displ
679
680	IF (field(ind)%ndim==2) THEN
681	ALLOCATE(Field_val2d(n))
682	n=0
683	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
684	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
685	k=d%iim*(j-1)+i
686	IF (d%own(i,j) .OR. single) THEN
687	n=n+1
688	Field_val2d(n)=field(ind)%rval2d(k)
689	ENDIF
690	ENDDO
691	ENDDO
692	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d, &
693	start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /))
694	DEALLOCATE(field_val2d)
695	ELSE IF (field(ind)%ndim==3) THEN
696	ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
697	n=0
698	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
699	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
700	k=d%iim*(j-1)+i
701	IF (d%own(i,j) .OR. single) THEN
702	n=n+1
703	Field_val3d(n,:)=field(ind)%rval3d(k,:)
704	ENDIF
705	ENDDO
706	ENDDO
707	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d, &
708	start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /))
709	DEALLOCATE(field_val3d)
710
711	ELSE IF (field(1)%ndim==4) THEN
712
713	DO q=1,FieldVarId(index)%size
714
715	ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
716	n=0
717	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
718	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
719	k=d%iim*(j-1)+i
720	IF (d%own(i,j) .OR. single) THEN
721	n=n+1
722	Field_val3d(n,:)=field(ind)%rval4d(k,:,q)
723	ENDIF
724	ENDDO
725	ENDDO
726
727	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d(:,l), &
728	start=(/ displ+ncell,l,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /))
729	DEALLOCATE(field_val3d)
730	ENDDO
731	ENDIF
732
733	ncell=ncell+n
734
735	ENDDO
736
737	ELSE IF (Field(ind_b)%field_type==field_Z) THEN
738	ncell=1
739	n=0
740	DO ind=ind_b,ind_e
741	d=>domain(ind)
742	CALL swap_geometry(ind)
743	CALL swap_dimensions(ind)
744
745	n=0
746	DO j=jj_begin+1,jj_end
747	DO i=ii_begin,ii_end-1
748	n=n+1
749	ENDDO
750	ENDDO
751
752	DO j=jj_begin,jj_end-1
753	DO i=ii_begin+1,ii_end
754	n=n+1
755	ENDDO
756	ENDDO
757
758	displ=FieldVarId(index)%displ
759
760	IF (field(ind)%ndim==2) THEN
761	ALLOCATE(Field_val2d(n))
762
763	n=0
764	DO j=jj_begin+1,jj_end
765	DO i=ii_begin,ii_end-1
766	n=n+1
767	k=iim*(j-1)+i
768	Field_val2d(n)=field(ind)%rval2d(k+z_down)
769	ENDDO
770	ENDDO
771
772	DO j=jj_begin,jj_end-1
773	DO i=ii_begin+1,ii_end
774	n=n+1
775	k=iim*(j-1)+i
776	Field_val2d(n)=field(ind)%rval2d(k+z_up)
777	ENDDO
778	ENDDO
779
780	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1), &
781	Field_val2d,start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /))
782	DEALLOCATE(field_val2d)
783
784	ELSE IF (field(ind)%ndim==3) THEN
785	ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
786	n=0
787	DO j=jj_begin+1,jj_end
788	DO i=ii_begin,ii_end-1
789	n=n+1
790	k=iim*(j-1)+i
791	Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:)
792	ENDDO
793	ENDDO
794
795	DO j=jj_begin,jj_end-1
796	DO i=ii_begin+1,ii_end
797	n=n+1
798	k=iim*(j-1)+i
799	Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:)
800	ENDDO
801	ENDDO
802
803	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d, &
804	start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /))
805	DEALLOCATE(field_val3d)
806
807	ELSE IF (field(1)%ndim==4) THEN
808
809	DO q=1,FieldVarId(index)%size
810	ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
811	n=0
812	DO j=jj_begin+1,jj_end
813	DO i=ii_begin,ii_end-1
814	n=n+1
815	k=iim*(j-1)+i
816	Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q)
817	ENDDO
818	ENDDO
819
820	DO j=jj_begin,jj_end-1
821	DO i=ii_begin+1,ii_end
822	n=n+1
823	k=iim*(j-1)+i
824	Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q)
825	ENDDO
826	ENDDO
827
828	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d, &
829	start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /))
830	DEALLOCATE(field_val3d)
831	ENDDO
832	ENDIF
833
834	ncell=ncell+n
835
836	ENDDO
837
838	ENDIF
839	status=NF90_SYNC(FieldId(Index))
840
841	END SUBROUTINE Writefield_mpi
842
843
844
845	! SUBROUTINE Create_header(name,field,nind)
846	! USE netcdf
847	! USE field_mod
848	! USE domain_mod
849	! USE spherical_geom_mod
850	! USE dimensions
851	! USE geometry
852	! IMPLICIT NONE
853	! CHARACTER(LEN=*) :: name
854	! TYPE(t_field),POINTER :: field(:)
855	! INTEGER,OPTIONAL,INTENT(IN) :: nind
856	! INTEGER :: ncell
857	! INTEGER :: nvert
858	! REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
859	! TYPE(t_domain),POINTER :: d
860	! INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
861	! INTEGER :: dim3id,dim4id
862	! INTEGER :: status
863	! INTEGER :: ind,i,j,k,n,q
864	! INTEGER :: iie,jje,iin,jjn
865	! INTEGER :: ind_b,ind_e
866	! INTEGER :: halo_size
867	! LOGICAL :: single
868	! INTEGER :: nij
869	!
870	! NbField=NbField+1
871	! FieldName(NbField)=TRIM(ADJUSTL(name))
872	! FieldIndex(NbField)=1
873	!
874	! IF (PRESENT(nind)) THEN
875	! ind_b=nind
876	! ind_e=nind
877	! halo_size=1
878	! single=.TRUE.
879	! ELSE
880	! ind_b=1
881	! ind_e=ndomain
882	! halo_size=0
883	! single=.FALSE.
884	! ENDIF
885	!
886	! ncell=0
887	!
888	! IF (Field(ind_b)%field_type==field_T) THEN
889	! nvert=6
890	!
891	! DO ind=ind_b,ind_e
892	! d=>domain(ind)
893	!
894	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
895	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
896	! IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1
897	! ENDDO
898	! ENDDO
899
900	! END DO
901	!
902	! status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
903	! FieldId(NbField)=ncid
904	! status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId)
905	! status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
906
907	! IF (Field(ind_b)%ndim==2) THEN
908	! FieldVarId(NbField)%size=1
909	! ALLOCATE(FieldVarId(NbField)%nc_id(1))
910	! ELSE IF (Field(ind_b)%ndim==3) THEN
911	! FieldVarId(NbField)%size=1
912	! ALLOCATE(FieldVarId(NbField)%nc_id(1))
913	! status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id)
914	! ELSE IF (Field(1)%ndim==4) THEN
915	! FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
916	! ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
917	! status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id)
918	! status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
919	! ENDIF
920	!
921	! status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
922	!
923	! status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
924	! status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
925	! status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
926	! status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
927	! status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
928	! status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
929	! status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
930	! status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
931	! status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
932	! status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
933
934	! IF (Field(ind_b)%ndim==2) THEN
935	! status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
936	! status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
937	! ELSE IF (Field(ind_b)%ndim==3) THEN
938	! status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
939	! status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
940	! ELSE IF (Field(ind_b)%ndim==4) THEN
941	! DO i=1,FieldVarId(NbField)%size
942	! status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),NF90_DOUBLE,(/ ncellId,dim3id,timeId /), &
943	! FieldVarId(NbField)%nc_id(i))
944	! status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat")
945	! ENDDO
946	! ENDIF
947	!
948	!
949	! status = NF90_ENDDEF(ncid)
950
951	! ncell=1
952	! DO ind=ind_b,ind_e
953	! d=>domain(ind)
954	!
955	! n=0
956	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
957	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
958	! IF (single .OR. d%own(i,j)) n=n+1
959	! ENDDO
960	! ENDDO
961	!
962	! ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
963	!
964	! n=0
965	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
966	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
967	! IF (d%own(i,j) .OR. single) THEN
968	! n=n+1
969	! CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n))
970	! lon(n)=lon(n)*180/Pi
971	! lat(n)=lat(n)*180/Pi
972	! DO k=0,5
973	! CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n))
974	! bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
975	! bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
976	! ENDDO
977	! ENDIF
978	! ENDDO
979	! ENDDO
980	! status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /))
981	! status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /))
982	! status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
983	! status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
984	!
985	! ncell=ncell+n
986	! DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
987	! END DO
988
989	! ELSE IF (Field(ind_b)%field_type==field_Z) THEN
990	! nvert=3
991	! DO ind=ind_b,ind_e
992	! d=>domain(ind)
993	!
994	! DO j=d%jj_begin+1,d%jj_end
995	! DO i=d%ii_begin,d%ii_end-1
996	! ncell=ncell+1
997	! ENDDO
998	! ENDDO
999
1000	! DO j=d%jj_begin,d%jj_end-1
1001	! DO i=d%ii_begin+1,d%ii_end
1002	! ncell=ncell+1
1003	! ENDDO
1004	! ENDDO
1005
1006	! END DO
1007	!
1008	! status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
1009	! FieldId(NbField)=ncid
1010	! status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId)
1011	! status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
1012
1013	! IF (Field(ind_b)%ndim==2) THEN
1014	! FieldVarId(NbField)%size=1
1015	! ALLOCATE(FieldVarId(NbField)%nc_id(1))
1016	! ELSE IF (Field(ind_b)%ndim==3) THEN
1017	! FieldVarId(NbField)%size=1
1018	! ALLOCATE(FieldVarId(NbField)%nc_id(1))
1019	! status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id)
1020	! ELSE IF (Field(1)%ndim==4) THEN
1021	! FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
1022	! ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
1023	! status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id)
1024	! ENDIF
1025
1026
1027	!
1028	! status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
1029	!
1030	! status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
1031	! status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
1032	! status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
1033	! status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
1034	! status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
1035	! status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
1036	! status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
1037	! status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
1038	! status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
1039	! status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
1040
1041
1042	! IF (Field(ind_b)%ndim==2) THEN
1043	! status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
1044	! status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
1045	! ELSE IF (Field(ind_b)%ndim==3) THEN
1046	! status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
1047	! status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
1048	! ELSE IF (Field(ind_b)%ndim==4) THEN
1049	! DO q=1,FieldVarId(NbField)%size
1050	! status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),NF90_DOUBLE, &
1051	! (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q))
1052	! status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat")
1053	! ENDDO
1054	! ENDIF
1055	!
1056	! status = NF90_ENDDEF(ncid)
1057
1058	! ncell=1
1059	! DO ind=ind_b,ind_e
1060	! d=>domain(ind)
1061	! CALL swap_geometry(ind)
1062	! CALL swap_dimensions(ind)
1063	!
1064	! n=0
1065	! DO j=jj_begin+1,jj_end
1066	! DO i=ii_begin,ii_end-1
1067	! n=n+1
1068	! ENDDO
1069	! ENDDO
1070
1071	! DO j=jj_begin,jj_end-1
1072	! DO i=ii_begin+1,ii_end
1073	! n=n+1
1074	! ENDDO
1075	! ENDDO
1076
1077	! ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
1078	!
1079	! n=0
1080	!
1081	! DO j=jj_begin+1,jj_end
1082	! DO i=ii_begin,ii_end-1
1083	! nij=(j-1)*iim+i
1084	! n=n+1
1085	! CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n))
1086	! lon(n)=lon(n)*180/Pi
1087	!! IF (lon(n)<0) lon(n)=lon(n)+360
1088	! lat(n)=lat(n)*180/Pi
1089	! CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
1090	! CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
1091	! CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
1092
1093	! DO k=0,2
1094	! bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1095	! bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1096	! IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
1097	! ENDDO
1098	! ENDDO
1099	! ENDDO
1100
1101	! DO j=jj_begin,jj_end-1
1102	! DO i=ii_begin+1,ii_end
1103	! nij=(j-1)*iim+i
1104	! n=n+1
1105	! CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
1106	! lon(n)=lon(n)*180/Pi
1107	! IF (lon(n)<0) lon(n)=lon(n)+360
1108	! lat(n)=lat(n)*180/Pi
1109	! CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
1110	! CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
1111	! CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
1112
1113	! DO k=0,2
1114	! bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1115	! bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1116	! IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
1117	! ENDDO
1118	! ENDDO
1119	! ENDDO
1120	!
1121	!
1122	! status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /))
1123	! status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /))
1124	! status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
1125	! status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
1126	!
1127	! ncell=ncell+n
1128	! DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
1129	! END DO
1130	! ENDIF
1131
1132
1133	!
1134	! status = NF90_CLOSE(ncid)
1135
1136	! END SUBROUTINE Create_Header
1137
1138
1139
1140
1141	SUBROUTINE Create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once)
1142	USE netcdf_mod
1143	USE field_mod
1144	USE domain_mod
1145	USE metric
1146	USE spherical_geom_mod
1147	IMPLICIT NONE
1148	CHARACTER(LEN=*),INTENT(IN) :: name_in
1149	TYPE(t_field),POINTER :: field(:)
1150	TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:)
1151	INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in
1152	INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in
1153	LOGICAL,INTENT(IN) :: once
1154
1155	INTEGER :: ncell
1156	INTEGER :: nvert
1157	REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
1158	TYPE(t_domain),POINTER :: d
1159	INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
1160	INTEGER :: dim3id,dim4id
1161	INTEGER :: status
1162	INTEGER :: ind,i,j,k,n,q
1163	INTEGER :: iie,jje,iin,jjn
1164	INTEGER :: ind_b,ind_e
1165	INTEGER :: halo_size
1166	LOGICAL :: single
1167	INTEGER :: nij
1168	CHARACTER(LEN=255) :: name
1169
1170
1171	name=TRIM(ADJUSTL(name_in))
1172
1173	IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN
1174	name=TRIM(name)//"_"//TRIM(int2str(ind_b))
1175	ind_b=ind_b_in
1176	ind_e=ind_b_in
1177	halo_size=1
1178	single=.TRUE.
1179	ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
1180	name=TRIM(name)//"_"//TRIM(int2str(ind_e))
1181	ind_b=ind_e_in
1182	ind_e=ind_e_in
1183	halo_size=1
1184	single=.TRUE.
1185	ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
1186	ind_b=ind_b_in
1187	ind_e=ind_e_in
1188	halo_size=0
1189	single=.FALSE.
1190	ELSE
1191	halo_size=0
1192	ind_b=1
1193	ind_e=ndomain
1194	single=.FALSE.
1195	ENDIF
1196
1197	NbField=NbField+1
1198	FieldName(NbField)=TRIM(ADJUSTL(name))
1199	FieldIndex(NbField)=1
1200
1201	ncell=0
1202
1203	IF (Field(ind_b)%field_type==field_T) THEN
1204	nvert=6
1205
1206	DO ind=ind_b,ind_e
1207	d=>domain_type(ind)
1208
1209	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
1210	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
1211	IF (single .OR. d%assign_domain(i,j)==ind) ncell=ncell+1
1212	ENDDO
1213	ENDDO
1214
1215	END DO
1216
1217	status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
1218	FieldId(NbField)=ncid
1219	status = NF90_DEF_DIM(ncid,'cell_i',ncell,ncellId)
1220	status = NF90_DEF_DIM(ncid,'nvert_i',nvert,nvertid)
1221
1222	IF (Field(ind_b)%ndim==2) THEN
1223	FieldVarId(NbField)%size=1
1224	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1225	ELSE IF (Field(ind_b)%ndim==3) THEN
1226	FieldVarId(NbField)%size=1
1227	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1228	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
1229	ELSE IF (Field(1)%ndim==4) THEN
1230	FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
1231	ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
1232	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
1233	! status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
1234	ENDIF
1235
1236	status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
1237
1238	status = NF90_DEF_VAR(ncid,'lon_i',NF90_DOUBLE,(/ ncellId /),lonId)
1239	status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
1240	status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
1241	status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon_i")
1242	status = NF90_DEF_VAR(ncid,'lat_i',NF90_DOUBLE,(/ ncellId /),latId)
1243	status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
1244	status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
1245	status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat_i")
1246	status = NF90_DEF_VAR(ncid,'bounds_lon_i',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
1247	status = NF90_DEF_VAR(ncid,'bounds_lat_i',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
1248
1249	IF (Field(ind_b)%ndim==2) THEN
1250	IF (once) THEN
1251	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId /),FieldVarId(NbField)%nc_id(1))
1252	ELSE
1253	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
1254	ENDIF
1255	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_i lat_i")
1256	ELSE IF (Field(ind_b)%ndim==3) THEN
1257	IF (once) THEN
1258	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(1))
1259	ELSE
1260	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
1261	ENDIF
1262	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_i lat_i")
1263	ELSE IF (Field(ind_b)%ndim==4) THEN
1264	DO i=1,FieldVarId(NbField)%size
1265	IF (once) THEN
1266	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id /), &
1267	FieldVarId(NbField)%nc_id(i))
1268	ELSE
1269	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /), &
1270	FieldVarId(NbField)%nc_id(i))
1271	ENDIF
1272	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon_i lat_i")
1273	ENDDO
1274	ENDIF
1275
1276
1277	status = NF90_ENDDEF(ncid)
1278
1279	ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell))
1280
1281	n=0
1282	DO ind=ind_b,ind_e
1283	d=>domain_type(ind)
1284	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
1285	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
1286	IF (d%assign_domain(i,j)==ind .OR. single) THEN
1287	n=n+1
1288	CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n))
1289	lon(n)=lon(n)*180/Pi
1290	lat(n)=lat(n)*180/Pi
1291	DO k=0,5
1292	CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n))
1293	bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1294	bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1295	ENDDO
1296	ENDIF
1297	ENDDO
1298	ENDDO
1299	ENDDO
1300	status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /))
1301	status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /))
1302	status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
1303	status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
1304
1305	DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
1306
1307	ELSE IF (Field(ind_b)%field_type==field_Z) THEN
1308	nvert=3
1309	DO ind=ind_b,ind_e
1310	d=>domain_type(ind)
1311
1312	DO j=d%jj_begin+1,d%jj_end
1313	DO i=d%ii_begin,d%ii_end-1
1314	ncell=ncell+1
1315	ENDDO
1316	ENDDO
1317
1318	DO j=d%jj_begin,d%jj_end-1
1319	DO i=d%ii_begin+1,d%ii_end
1320	ncell=ncell+1
1321	ENDDO
1322	ENDDO
1323
1324	END DO
1325
1326	status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
1327	FieldId(NbField)=ncid
1328	status = NF90_DEF_DIM(ncid,'cell_v',ncell,ncellId)
1329	status = NF90_DEF_DIM(ncid,'nvert_v',nvert,nvertid)
1330
1331	IF (Field(ind_b)%ndim==2) THEN
1332	FieldVarId(NbField)%size=1
1333	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1334	ELSE IF (Field(ind_b)%ndim==3) THEN
1335	FieldVarId(NbField)%size=1
1336	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1337	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
1338	ELSE IF (Field(1)%ndim==4) THEN
1339	FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
1340	ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
1341	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
1342	ENDIF
1343
1344
1345
1346	status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
1347
1348	status = NF90_DEF_VAR(ncid,'lon_v',NF90_DOUBLE,(/ ncellId /),lonId)
1349	status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
1350	status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
1351	status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon_v")
1352	status = NF90_DEF_VAR(ncid,'lat_v',NF90_DOUBLE,(/ ncellId /),latId)
1353	status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
1354	status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
1355	status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat_v")
1356	status = NF90_DEF_VAR(ncid,'bounds_lon_v',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
1357	status = NF90_DEF_VAR(ncid,'bounds_lat_v',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
1358
1359
1360	IF (Field(ind_b)%ndim==2) THEN
1361	IF (once) THEN
1362	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId /),FieldVarId(NbField)%nc_id(1))
1363	ELSE
1364	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
1365	ENDIF
1366	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_v lat_v")
1367	ELSE IF (Field(ind_b)%ndim==3) THEN
1368	IF (once) THEN
1369	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(1))
1370	ELSE
1371	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
1372	ENDIF
1373	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_v lat_v")
1374	ELSE IF (Field(ind_b)%ndim==4) THEN
1375	DO q=1,FieldVarId(NbField)%size
1376	IF (once) THEN
1377	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),ncprec, &
1378	(/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(q))
1379	ELSE
1380	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),ncprec, &
1381	(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q))
1382	ENDIF
1383	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon_v lat_v")
1384	ENDDO
1385	ENDIF
1386
1387	status = NF90_ENDDEF(ncid)
1388
1389	ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell))
1390
1391	n=0
1392
1393	DO ind=ind_b,ind_e
1394	d=>domain_type(ind)
1395	DO j=d%jj_begin+1,d%jj_end
1396	DO i=d%ii_begin,d%ii_end-1
1397	nij=(j-1)*d%iim+i
1398	n=n+1
1399	CALL xyz2lonlat(d%vertex(:,vdown,i,j),lon(n),lat(n))
1400	lon(n)=lon(n)*180/Pi
1401	lat(n)=lat(n)*180/Pi
1402
1403	CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n))
1404	CALL xyz2lonlat(d%xyz(:,i,j-1),bounds_lon(1,n), bounds_lat(1,n))
1405	CALL xyz2lonlat(d%xyz(:,i+1,j-1),bounds_lon(2,n), bounds_lat(2,n))
1406
1407	DO k=0,2
1408	bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1409	bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1410	ENDDO
1411	ENDDO
1412	ENDDO
1413
1414	DO j=d%jj_begin,d%jj_end-1
1415	DO i=d%ii_begin+1,d%ii_end
1416	nij=(j-1)*d%iim+i
1417	n=n+1
1418	CALL xyz2lonlat(d%vertex(:,vup,i,j),lon(n),lat(n))
1419	lon(n)=lon(n)*180/Pi
1420	lat(n)=lat(n)*180/Pi
1421	CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n))
1422	CALL xyz2lonlat(d%xyz(:,i,j+1),bounds_lon(1,n), bounds_lat(1,n))
1423	CALL xyz2lonlat(d%xyz(:,i-1,j+1),bounds_lon(2,n), bounds_lat(2,n))
1424
1425	DO k=0,2
1426	bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1427	bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1428	ENDDO
1429	ENDDO
1430	ENDDO
1431	ENDDO
1432
1433	status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /))
1434	status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /))
1435	status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
1436	status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
1437
1438	DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
1439	ENDIF
1440
1441
1442	END SUBROUTINE Create_Header_gen
1443
1444	SUBROUTINE Create_header_mpi(name,field,nind)
1445	USE netcdf_mod
1446	USE field_mod
1447	USE domain_mod
1448	USE spherical_geom_mod
1449	USE dimensions
1450	USE geometry
1451	USE mpi_mod
1452	USE mpipara
1453	IMPLICIT NONE
1454	CHARACTER(LEN=*) :: name
1455	CHARACTER(LEN=LEN_TRIM(ADJUSTL(name))) :: name_adj
1456	TYPE(t_field),POINTER :: field(:)
1457	INTEGER,OPTIONAL,INTENT(IN) :: nind
1458	INTEGER :: ncell
1459	INTEGER :: nvert
1460	REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
1461	TYPE(t_domain),POINTER :: d
1462	INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
1463	INTEGER :: dim3id,dim4id
1464	INTEGER :: status
1465	INTEGER :: ind,i,j,k,n,q
1466	INTEGER :: iie,jje,iin,jjn
1467	INTEGER :: ind_b,ind_e
1468	INTEGER :: halo_size
1469	LOGICAL :: single
1470	INTEGER :: nij
1471	INTEGER :: ncell_glo(0:mpi_size-1)
1472	INTEGER :: displ, ncell_tot
1473
1474
1475	NbField=NbField+1
1476	name_adj=TRIM(ADJUSTL(name)) ! work around ICE with pgf90
1477	FieldName(NbField)=name_adj
1478	FieldIndex(NbField)=1
1479
1480	IF (PRESENT(nind)) THEN
1481	ind_b=nind
1482	ind_e=nind
1483	halo_size=1
1484	single=.TRUE.
1485	ELSE
1486	ind_b=1
1487	ind_e=ndomain
1488	halo_size=0
1489	single=.FALSE.
1490	ENDIF
1491
1492	ncell=0
1493
1494	IF (Field(ind_b)%field_type==field_T) THEN
1495	nvert=6
1496
1497	DO ind=ind_b,ind_e
1498	d=>domain(ind)
1499
1500	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
1501	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
1502	IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1
1503	ENDDO
1504	ENDDO
1505
1506	END DO
1507
1508	CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr)
1509
1510	displ=0
1511	DO i=1,mpi_rank
1512	displ=displ+ncell_glo(i-1)
1513	ENDDO
1514	FieldVarId(NbField)%displ=displ
1515	ncell_tot=sum(ncell_glo(:))
1516
1517	! status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc', IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid)
1518	FieldId(NbField)=ncid
1519	status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId)
1520	status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
1521
1522	IF (Field(ind_b)%ndim==2) THEN
1523	FieldVarId(NbField)%size=1
1524	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1525	ELSE IF (Field(ind_b)%ndim==3) THEN
1526	FieldVarId(NbField)%size=1
1527	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1528	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
1529	ELSE IF (Field(1)%ndim==4) THEN
1530	FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
1531	ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
1532	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
1533	! status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
1534	ENDIF
1535
1536	status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
1537
1538	status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
1539	status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
1540	status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
1541	status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
1542	status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
1543	status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
1544	status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
1545	status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
1546	status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
1547	status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
1548
1549	IF (Field(ind_b)%ndim==2) THEN
1550	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
1551	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
1552	status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/))
1553	ELSE IF (Field(ind_b)%ndim==3) THEN
1554	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
1555	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
1556	status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, &
1557	(/ncell_tot,size(field(ind_b)%rval3d,2),1/))
1558	ELSE IF (Field(ind_b)%ndim==4) THEN
1559	DO i=1,FieldVarId(NbField)%size
1560	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /), &
1561	FieldVarId(NbField)%nc_id(i))
1562	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat")
1563	status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED, &
1564	(/ncell_tot,size(field(ind_b)%rval4d,2),1/))
1565	ENDDO
1566	ENDIF
1567
1568
1569	status = NF90_ENDDEF(ncid)
1570
1571	ncell=1
1572	DO ind=ind_b,ind_e
1573	d=>domain(ind)
1574
1575	n=0
1576	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
1577	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
1578	IF (single .OR. d%own(i,j)) n=n+1
1579	ENDDO
1580	ENDDO
1581
1582	ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
1583
1584	n=0
1585	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
1586	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
1587	IF (d%own(i,j) .OR. single) THEN
1588	n=n+1
1589	CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n))
1590	lon(n)=lon(n)*180/Pi
1591	lat(n)=lat(n)*180/Pi
1592	DO k=0,5
1593	CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n))
1594	bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1595	bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1596	ENDDO
1597	ENDIF
1598	ENDDO
1599	ENDDO
1600	status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /))
1601	status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /))
1602	status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
1603	status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
1604
1605	ncell=ncell+n
1606	DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
1607	END DO
1608
1609	ELSE IF (Field(ind_b)%field_type==field_Z) THEN
1610	nvert=3
1611	DO ind=ind_b,ind_e
1612	d=>domain(ind)
1613
1614	DO j=d%jj_begin+1,d%jj_end
1615	DO i=d%ii_begin,d%ii_end-1
1616	ncell=ncell+1
1617	ENDDO
1618	ENDDO
1619
1620	DO j=d%jj_begin,d%jj_end-1
1621	DO i=d%ii_begin+1,d%ii_end
1622	ncell=ncell+1
1623	ENDDO
1624	ENDDO
1625
1626	END DO
1627
1628	CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr)
1629
1630	displ=0
1631	DO i=1,mpi_rank
1632	displ=displ+ncell_glo(i-1)
1633	ENDDO
1634	FieldVarId(NbField)%displ=displ
1635	ncell_tot=sum(ncell_glo(:))
1636
1637	! status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc',IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid)
1638	! status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
1639	FieldId(NbField)=ncid
1640	status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId)
1641	status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
1642
1643	IF (Field(ind_b)%ndim==2) THEN
1644	FieldVarId(NbField)%size=1
1645	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1646	ELSE IF (Field(ind_b)%ndim==3) THEN
1647	FieldVarId(NbField)%size=1
1648	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1649	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
1650	ELSE IF (Field(1)%ndim==4) THEN
1651	FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
1652	ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
1653	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
1654	ENDIF
1655
1656
1657
1658	status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
1659
1660	status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
1661	status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
1662	status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
1663	status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
1664	status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
1665	status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
1666	status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
1667	status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
1668	status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
1669	status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
1670
1671
1672	IF (Field(ind_b)%ndim==2) THEN
1673	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
1674	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
1675	status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/))
1676	ELSE IF (Field(ind_b)%ndim==3) THEN
1677	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
1678	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
1679	status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, &
1680	(/ncell_tot,size(field(ind_b)%rval3d,2),1/))
1681	ELSE IF (Field(ind_b)%ndim==4) THEN
1682	DO q=1,FieldVarId(NbField)%size
1683	status = NF90_DEF_VAR(ncid,name_adj//int2str(q),ncprec, &
1684	(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q))
1685	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat")
1686	status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED, &
1687	(/ncell_tot,size(field(ind_b)%rval4d,2),1/))
1688	ENDDO
1689	ENDIF
1690
1691	status = NF90_ENDDEF(ncid)
1692
1693	ncell=1
1694	DO ind=ind_b,ind_e
1695	d=>domain(ind)
1696	CALL swap_geometry(ind)
1697	CALL swap_dimensions(ind)
1698
1699	n=0
1700	DO j=jj_begin+1,jj_end
1701	DO i=ii_begin,ii_end-1
1702	n=n+1
1703	ENDDO
1704	ENDDO
1705
1706	DO j=jj_begin,jj_end-1
1707	DO i=ii_begin+1,ii_end
1708	n=n+1
1709	ENDDO
1710	ENDDO
1711
1712	ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
1713
1714	n=0
1715
1716	DO j=jj_begin+1,jj_end
1717	DO i=ii_begin,ii_end-1
1718	nij=(j-1)*iim+i
1719	n=n+1
1720	CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n))
1721	lon(n)=lon(n)*180/Pi
1722	lat(n)=lat(n)*180/Pi
1723	CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
1724	CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
1725	CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
1726
1727	DO k=0,2
1728	bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1729	bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1730	ENDDO
1731	ENDDO
1732	ENDDO
1733
1734	DO j=jj_begin,jj_end-1
1735	DO i=ii_begin+1,ii_end
1736	nij=(j-1)*iim+i
1737	n=n+1
1738	CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
1739	lon(n)=lon(n)*180/Pi
1740	lat(n)=lat(n)*180/Pi
1741	CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
1742	CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
1743	CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
1744
1745	DO k=0,2
1746	bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1747	bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1748	ENDDO
1749	ENDDO
1750	ENDDO
1751
1752
1753	status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /))
1754	status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /))
1755	status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
1756	status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
1757
1758	ncell=ncell+n
1759	DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
1760	END DO
1761	ENDIF
1762
1763
1764	END SUBROUTINE Create_Header_mpi
1765
1766	SUBROUTINE Close_files
1767	USE netcdf
1768	IMPLICIT NONE
1769	INTEGER :: i,k,status
1770
1771	DO i=1,NbField
1772	status=NF90_CLOSE(FieldId(i))
1773	ENDDO
1774	END SUBROUTINE Close_files
1775
1776
1777	function int2str(int)
1778	implicit none
1779	integer, parameter :: MaxLen=10
1780	integer,intent(in) :: int
1781	character(len=MaxLen) :: int2str
1782	logical :: flag
1783	integer :: i
1784	flag=.true.
1785
1786	i=int
1787
1788	int2str=''
1789	do while (flag)
1790	int2str=CHAR(MOD(i,10)+48)//int2str
1791	i=i/10
1792	if (i==0) flag=.false.
1793	enddo
1794	end function int2str
1795
1796	end module write_field
1797

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

Download in other formats: