Context Navigation

source: codes/icosagcm/trunk/src/output/write_field.f90 @ 581

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

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

Download in other formats: