Context Navigation

write_field.f90 @ 3891

Last change on this file since 3891 was 3861, checked in by ymipsl, 9 years ago

Replace XIOS1 by XIOS2 in DYNAMICO

File size: 63.3 KB

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

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

Download in other formats:

Original Format