Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

write_field.f90

Last change on this file was 963, checked in by adurocher, 5 years ago
Merge 'mpi_rewrite' into trunk
File size: 63.1 KB

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

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

Context Navigation

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

Download in other formats: