Context Navigation

source: codes/icosagcm/devel/src/output/write_field.f90 @ 879

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

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

Download in other formats: