Context Navigation

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

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

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

Download in other formats: