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

Last change on this file since 880 was 880, checked in by dubos, 5 years ago

devel : store cell bounds once, use them for XIOS later
File size: 15.6 KB
Line 
1module write_field_mod
2  USE genmod
3  USE write_field_vars_mod
4  IMPLICIT NONE
5  PRIVATE 
6 
7  PUBLIC init_writeField, writefield, close_files
8 
9  CONTAINS
10 
11    SUBROUTINE init_writeField
12      USE ioipsl
13      use netcdf_mod
14      CHARACTER(LEN=255) :: netcdf_prec
15     
16      !$OMP CRITICAL
17      netcdf_prec='float'
18      CALL getin("netcdf_prec",netcdf_prec)
19      SELECT CASE(TRIM(netcdf_prec))
20        CASE('float')
21          ncprec=NF90_FLOAT
22        CASE('double')
23          ncprec=NF90_DOUBLE
24        CASE default
25        PRINT*,'Bad selector for variable netcdf_prec : <', TRIM(netcdf_prec),"> options are <float>, <double>" 
26        STOP   
27      END SELECT
28      !$OMP END CRITICAL
29    END SUBROUTINE init_writeField
30   
31   SUBROUTINE Writefield(name_in,field,nind,once)
32    USE domain_mod
33    USE field_mod
34    USE transfert_mpi_mod
35    USE dimensions
36    USE mpipara
37    USE netcdf_mod
38    USE grid_param
39     CHARACTER(LEN=*),INTENT(IN) :: name_in
40      TYPE(t_field),POINTER :: field(:)
41      INTEGER,OPTIONAL,INTENT(IN) :: nind
42      LOGICAL,OPTIONAL,INTENT(IN) :: once
43      LOGICAL                     :: once_
44     
45      TYPE(t_field),POINTER :: field_glo(:)
46     
47      IF (no_io) RETURN
48
49!$OMP BARRIER
50!$OMP MASTER     
51      IF(PRESENT(once)) THEN
52         once_=once
53      ELSE
54         once_=.FALSE.
55      END IF
56
57      IF (grid_type==grid_ico) THEN
58         IF (field(1)%ndim==2) THEN
59            CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type)
60         ELSE IF (field(1)%ndim==3) THEN
61            CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3) 
62         ELSE IF (field(1)%ndim==4) THEN
63            CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3,field(1)%dim4)
64         ENDIF
65         
66         CALL gather_field(field,field_glo)
67      ELSE
68         field_glo => field ! FIXME unstructured + MPI
69      END IF
70
71      IF (mpi_rank==0) THEN
72        IF (PRESENT(nind)) THEN
73         CALL writefield_gen(name_in,field_glo,domain_glo,nind,once=once_)
74        ELSE
75         CALL writefield_gen(name_in,field_glo,domain_glo,1,ndomain_glo,once=once_)
76        ENDIF
77      ENDIF
78     
79      IF(grid_type == grid_ico) CALL deallocate_field_glo(field_glo)
80!$OMP END MASTER
81!$OMP BARRIER
82     
83   END SUBROUTINE writefield
84   
85   SUBROUTINE Writefield_gen(name_in, field, domain_type, ind_b_in, ind_e_in,once )
86     USE netcdf_mod
87     USE domain_mod
88     USE field_mod
89     CHARACTER(LEN=*),INTENT(IN) :: name_in
90     TYPE(t_field), POINTER :: field(:)
91     TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:)
92     INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in
93     INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in
94     LOGICAL, INTENT(IN) :: once
95     ! local variables
96     TYPE(t_cellset), POINTER :: cells(:)
97     TYPE(t_domain),POINTER :: d
98     LOGICAL :: single
99     INTEGER :: ind_b, ind_e, ind
100     CHARACTER(len=255) :: name
101     INTEGER :: index, ncell, nvert, n
102     ! for embedded routines
103     REAL(r8),ALLOCATABLE :: field_val2d(:)
104     REAL(r8),ALLOCATABLE :: field_val3d(:,:)
105     INTEGER :: status, n_begin, ij, q, dim3
106
107     name=TRIM(ADJUSTL(name_in))
108
109      IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN
110        name=TRIM(name)//"_"//TRIM(int2str(ind_b))
111        ind_b=ind_b_in
112        ind_e=ind_b_in
113        single=.TRUE.
114      ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
115        name=TRIM(name)//"_"//TRIM(int2str(ind_e))
116        ind_b=ind_e_in
117        ind_e=ind_e_in
118        single=.TRUE.
119      ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
120        ind_b=ind_b_in
121        ind_e=ind_e_in
122        single=.FALSE.
123      ELSE
124        ind_b=1
125        ind_e=ndomain
126        single=.FALSE.
127      ENDIF     
128
129      index=GetFieldIndex(name)
130      if (index==-1) then
131         call create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once)
132         index=GetFieldIndex(name)
133      else
134         FieldIndex(index)=FieldIndex(index)+1.
135      endif
136
137      SELECT CASE(Field(ind_b)%field_type)
138      CASE(field_T)
139         IF(single) THEN ! include halos
140            cells => mesh_glo%primal_all
141         ELSE
142            cells => mesh_glo%primal_own
143         END IF
144      CASE(field_Z)
145         IF(single) THEN ! include halos
146            cells => mesh_glo%dual_all
147         ELSE
148            cells => mesh_glo%dual_own
149         END IF
150      END SELECT
151
152      ncell=0
153      DO ind=ind_b,ind_e
154         nvert=SIZE(cells(ind)%bnds_lon,1)
155         ncell = ncell + cells(ind)%ncell
156      END DO
157
158      SELECT CASE(field(1)%ndim)
159      CASE(2)
160         CALL write_2d
161      CASE(3)
162         CALL write_3d
163      CASE(4)         
164         CALL write_4d
165      END SELECT
166     
167      status=NF90_SYNC(FieldId(index))
168     
169   CONTAINS
170   
171     SUBROUTINE write_2d
172       ALLOCATE(Field_val2d(ncell))
173       n_begin=0 
174       DO ind=ind_b,ind_e
175          DO n=1, cells(ind)%ncell
176             ij=cells(ind)%ij(n)
177             PRINT *, 'write_2d :', ind, n, n_begin+n, ij 
178             field_val2d(n_begin+n) = field(ind)%rval2d(ij)
179          END DO
180          n_begin = n_begin + cells(ind)%ncell
181       END DO
182       IF (once) THEN
183          status=NF90_PUT_VAR(FieldId(index), FieldVarId(index)%nc_id(1), &
184               Field_val2d, start=(/ 1 /),count=(/ncell /))
185       ELSE
186          status=NF90_PUT_VAR(FieldId(Index), FieldVarId(index)%nc_id(1), &
187               Field_val2d,  start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /))
188       ENDIF
189       DEALLOCATE(field_val2d)
190     END SUBROUTINE write_2d
191
192     SUBROUTINE write_3d
193       dim3 = SIZE(field(1)%rval3d,2)
194       ALLOCATE(field_val3d(ncell,dim3))
195       n_begin=0 
196       DO ind=ind_b,ind_e
197          DO n=1, cells(ind)%ncell
198             ij=cells(ind)%ij(n)
199             field_val3d(n_begin+n,:) = field(ind)%rval3d(ij,:)
200          END DO
201       END DO
202       IF (once) THEN
203          status=NF90_PUT_VAR(FieldId(Index), FieldVarId(index)%nc_id(1), &
204               field_val3d, start=(/ 1,1 /), &
205               count=(/ncell,dim3 /))
206       ELSE
207          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1), &
208               field_val3d, start=(/ 1,1,FieldIndex(Index) /),   &
209               count=(/ncell,size(field(1)%rval3d,2),1 /))
210       ENDIF
211       DEALLOCATE(field_val3d)
212     END SUBROUTINE write_3d
213
214     SUBROUTINE write_4d
215       dim3 = SIZE(field(1)%rval4d,2)
216       ALLOCATE(field_val3d(ncell,dim3))
217       DO q=1,FieldVarId(index)%size
218          n_begin=0 
219          DO ind=ind_b,ind_e
220             DO n=1, cells(ind)%ncell
221                ij=cells(ind)%ij(n)
222                field_val3d(n_begin+n,:) = field(ind)%rval4d(ij,:,q)
223             END DO
224          END DO
225          IF (once) THEN
226             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q), &
227                  field_val3d, start=(/ 1,1 /),   &
228                  count=(/ncell,dim3,1 /))
229          ELSE
230             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q), &
231                  field_val3d, start=(/ 1,1,FieldIndex(Index) /),   &
232                  count=(/ncell,dim3,1 /))
233          ENDIF
234       END DO
235       DEALLOCATE(field_val3d)
236     END SUBROUTINE write_4d
237
238   END SUBROUTINE Writefield_gen
239   
240   SUBROUTINE Create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once)
241     USE netcdf_mod
242     USE field_mod
243     USE domain_mod
244     USE metric
245     USE spherical_geom_mod
246     CHARACTER(LEN=*),INTENT(IN) :: name_in
247     TYPE(t_field),POINTER :: field(:)
248     TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:)
249     INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in
250     INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in
251     LOGICAL,INTENT(IN) :: once
252     
253     TYPE(t_cellset), POINTER :: cells(:)
254     INTEGER :: ncell
255     INTEGER :: nvert
256     REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
257     TYPE(t_domain),POINTER :: d
258     INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
259     INTEGER :: dim3id,dim4id
260     INTEGER :: status
261     INTEGER :: ind,i,j,k,n,q
262     INTEGER :: iie,jje,iin,jjn
263     INTEGER :: ind_b,ind_e, n_begin, n_end
264     LOGICAL :: single 
265     INTEGER :: nij
266     CHARACTER(LEN=255) :: name
267     INTEGER :: l,level_size, levId, dimlevId
268     
269     name=TRIM(ADJUSTL(name_in))
270     
271     IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN
272        name=TRIM(name)//"_"//TRIM(int2str(ind_b))
273        ind_b=ind_b_in
274        ind_e=ind_b_in
275        single=.TRUE.
276     ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
277        name=TRIM(name)//"_"//TRIM(int2str(ind_e))
278        ind_b=ind_e_in
279        ind_e=ind_e_in
280        single=.TRUE.
281     ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
282        ind_b=ind_b_in
283        ind_e=ind_e_in
284        single=.FALSE.
285     ELSE
286        ind_b=1
287        ind_e=ndomain
288        single=.FALSE.
289     ENDIF
290     
291     NbField=NbField+1
292     FieldName(NbField)=TRIM(ADJUSTL(name))
293     FieldIndex(NbField)=1
294     
295     PRINT *, 'Creating header for writefield : ', TRIM(FieldName(NbField)), ndomain, ind_b, ind_e ! FIXME
296     PRINT *, mesh_glo%ndomain, SIZE(mesh_glo%primal_own)
297     SELECT CASE(Field(ind_b)%field_type)
298     CASE(field_T)
299        PRINT *, '    field_type == field_T' ! FIXME
300        IF(single) THEN ! include halos
301           cells => mesh_glo%primal_all
302        ELSE
303           cells => mesh_glo%primal_own
304        END IF
305     CASE(field_Z)
306        PRINT *, '    field_type == field_Z' ! FIXME
307        IF(single) THEN ! include halos
308           cells => mesh_glo%dual_all
309        ELSE
310           cells => mesh_glo%dual_own
311        END IF
312     END SELECT
313     
314     PRINT *, 'writefield : ind_b, ind_e :', ind_b, ind_e
315
316     ncell=0
317     DO ind=ind_b,ind_e
318        nvert=SIZE(cells(ind)%bnds_lon,1)
319        ncell = ncell + cells(ind)%ncell
320     END DO
321     PRINT *, 'writefield : found ',ncell, ' cells.' ! FIXME
322     
323     status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
324     FieldId(NbField)=ncid
325     status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
326     
327     level_size=0
328     SELECT CASE(Field(ind_b)%ndim)
329     CASE(2)
330        FieldVarId(NbField)%size=1
331        ALLOCATE(FieldVarId(NbField)%nc_id(1))
332     CASE(3)
333        FieldVarId(NbField)%size=1
334        ALLOCATE(FieldVarId(NbField)%nc_id(1))
335        status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
336        level_size=size(field(ind_b)%rval3d,2)
337     CASE(4)
338        FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
339        ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
340        status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
341        level_size=size(field(ind_b)%rval4d,2)
342        !          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
343     END SELECT
344     
345     PRINT*,"create_header_gen : LEVEL_SIZE=",level_size
346     IF (level_size>0) THEN
347        status = NF90_DEF_VAR(ncid,'lev',NF90_DOUBLE,(/ dim3id /),levId)
348        status = NF90_PUT_ATT(ncid,levId,"axis","Z")
349     ENDIF
350     
351     SELECT CASE(Field(ind_b)%field_type)
352     CASE(field_T)
353        status = NF90_DEF_DIM(ncid,'cell_i',ncell,ncellId)
354        status = NF90_DEF_DIM(ncid,'nvert_i',nvert,nvertid)
355        status = NF90_DEF_VAR(ncid,'lon_i',NF90_DOUBLE,(/ ncellId /),lonId)
356        status = NF90_DEF_VAR(ncid,'lat_i',NF90_DOUBLE,(/ ncellId /),latId)
357        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon_i")
358        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat_i")
359        status = NF90_DEF_VAR(ncid,'bounds_lon_i',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
360        status = NF90_DEF_VAR(ncid,'bounds_lat_i',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
361     CASE(field_Z)
362        status = NF90_DEF_DIM(ncid,'cell_v',ncell,ncellId)
363        status = NF90_DEF_DIM(ncid,'nvert_v',nvert,nvertid)
364        status = NF90_DEF_VAR(ncid,'lon_v',NF90_DOUBLE,(/ ncellId /),lonId)
365        status = NF90_DEF_VAR(ncid,'lat_v',NF90_DOUBLE,(/ ncellId /),latId)
366        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon_v")
367        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat_v")
368        status = NF90_DEF_VAR(ncid,'bounds_lon_v',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
369        status = NF90_DEF_VAR(ncid,'bounds_lat_v',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
370     END SELECT
371     
372     status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
373     status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
374     status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
375     status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
376     
377     SELECT CASE(Field(ind_b)%ndim)
378     CASE(2)
379        IF (once) THEN
380           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, &
381                (/ ncellId /),FieldVarId(NbField)%nc_id(1))
382        ELSE
383           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, &
384                (/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
385        END IF
386        CALL put_att_coordinates(1)
387     CASE(3)
388        IF (once) THEN
389           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, &
390                (/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(1))
391        ELSE
392           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, &
393                (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
394        END IF
395        CALL put_att_coordinates(1)
396     CASE(4)
397        DO i=1,FieldVarId(NbField)%size
398           IF (once) THEN
399              status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id /),  &
400                   FieldVarId(NbField)%nc_id(i))
401           ELSE
402              status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /),  &
403                   FieldVarId(NbField)%nc_id(i))
404           END IF
405           CALL put_att_coordinates(i)
406           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon_i lat_i")
407        END DO
408     END SELECT
409     
410     status = NF90_ENDDEF(ncid)     
411     
412     if (level_size>0) status = NF90_PUT_VAR(ncid,levId,(/ (l,l=1,level_size) /))
413     
414     ALLOCATE(lon(ncell),lat(ncell))
415     ALLOCATE(bounds_lon(nvert,ncell), bounds_lat(nvert,ncell))
416     n_begin=0 
417     DO ind=ind_b,ind_e
418        n_end = n_begin + cells(ind)%ncell
419        PRINT *, 'create_header_gen ', n_begin, n_end, SHAPE(cells(ind)%bnds_lon), SHAPE(bounds_lon)
420        lat(n_begin+1:n_end) = cells(ind)%lat(:) *180./Pi
421        lon(n_begin+1:n_end) = cells(ind)%lon(:) *180./Pi
422        bounds_lon(:,n_begin+1:n_end) = cells(ind)%bnds_lon(:,:) *180./Pi
423        bounds_lat(:,n_begin+1:n_end) = cells(ind)%bnds_lat(:,:) *180./Pi
424        n_begin = n_end
425     END DO
426     status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /))
427     status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /))
428     status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
429     status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
430     
431     DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
432     
433   CONTAINS
434     
435     SUBROUTINE put_att_coordinates(ind)
436       INTEGER :: ind
437       SELECT CASE(Field(ind_b)%field_type)
438       CASE(field_T)
439          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(ind), &
440               "coordinates","lon_i lat_i")
441       CASE(field_Z)
442          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(ind), &
443               "coordinates","lon_v lat_v")     
444       END SELECT
445     END SUBROUTINE put_att_coordinates
446     
447   END SUBROUTINE Create_Header_gen
448   
449   SUBROUTINE Close_files
450   USE netcdf_mod
451     INTEGER :: i,k,status
452!$OMP MASTER     
453     DO i=1,NbField
454         status=NF90_CLOSE(FieldId(i))
455     ENDDO
456!$OMP END MASTER
457   END SUBROUTINE  Close_files
458     
459     
460end module write_field_mod
Note: See TracBrowser for help on using the repository browser.