Context Navigation

write_field.f90 @ 3845

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

Dynamico side modification for Dynamico-LMDZ coupling

File size: 62.8 KB

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

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

Download in other formats:

Original Format