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

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

devel : quiet down write_field.f90
File size: 15.5 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             field_val2d(n_begin+n) = field(ind)%rval2d(ij)
178          END DO
179          n_begin = n_begin + cells(ind)%ncell
180       END DO
181       IF (once) THEN
182          status=NF90_PUT_VAR(FieldId(index), FieldVarId(index)%nc_id(1), &
183               Field_val2d, start=(/ 1 /),count=(/ncell /))
184       ELSE
185          status=NF90_PUT_VAR(FieldId(Index), FieldVarId(index)%nc_id(1), &
186               Field_val2d,  start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /))
187       ENDIF
188       DEALLOCATE(field_val2d)
189     END SUBROUTINE write_2d
190
191     SUBROUTINE write_3d
192       dim3 = SIZE(field(1)%rval3d,2)
193       ALLOCATE(field_val3d(ncell,dim3))
194       n_begin=0 
195       DO ind=ind_b,ind_e
196          DO n=1, cells(ind)%ncell
197             ij=cells(ind)%ij(n)
198             field_val3d(n_begin+n,:) = field(ind)%rval3d(ij,:)
199          END DO
200       END DO
201       IF (once) THEN
202          status=NF90_PUT_VAR(FieldId(Index), FieldVarId(index)%nc_id(1), &
203               field_val3d, start=(/ 1,1 /), &
204               count=(/ncell,dim3 /))
205       ELSE
206          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1), &
207               field_val3d, start=(/ 1,1,FieldIndex(Index) /),   &
208               count=(/ncell,size(field(1)%rval3d,2),1 /))
209       ENDIF
210       DEALLOCATE(field_val3d)
211     END SUBROUTINE write_3d
212
213     SUBROUTINE write_4d
214       dim3 = SIZE(field(1)%rval4d,2)
215       ALLOCATE(field_val3d(ncell,dim3))
216       DO q=1,FieldVarId(index)%size
217          n_begin=0 
218          DO ind=ind_b,ind_e
219             DO n=1, cells(ind)%ncell
220                ij=cells(ind)%ij(n)
221                field_val3d(n_begin+n,:) = field(ind)%rval4d(ij,:,q)
222             END DO
223          END DO
224          IF (once) THEN
225             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q), &
226                  field_val3d, start=(/ 1,1 /),   &
227                  count=(/ncell,dim3,1 /))
228          ELSE
229             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q), &
230                  field_val3d, start=(/ 1,1,FieldIndex(Index) /),   &
231                  count=(/ncell,dim3,1 /))
232          ENDIF
233       END DO
234       DEALLOCATE(field_val3d)
235     END SUBROUTINE write_4d
236
237   END SUBROUTINE Writefield_gen
238   
239   SUBROUTINE Create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once)
240     USE netcdf_mod
241     USE field_mod
242     USE domain_mod
243     USE metric
244     USE spherical_geom_mod
245     CHARACTER(LEN=*),INTENT(IN) :: name_in
246     TYPE(t_field),POINTER :: field(:)
247     TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:)
248     INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in
249     INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in
250     LOGICAL,INTENT(IN) :: once
251     
252     TYPE(t_cellset), POINTER :: cells(:)
253     INTEGER :: ncell
254     INTEGER :: nvert
255     REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
256     TYPE(t_domain),POINTER :: d
257     INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
258     INTEGER :: dim3id,dim4id
259     INTEGER :: status
260     INTEGER :: ind,i,j,k,n,q
261     INTEGER :: iie,jje,iin,jjn
262     INTEGER :: ind_b,ind_e, n_begin, n_end
263     LOGICAL :: single 
264     INTEGER :: nij
265     CHARACTER(LEN=255) :: name
266     INTEGER :: l,level_size, levId, dimlevId
267     
268     name=TRIM(ADJUSTL(name_in))
269     
270     IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN
271        name=TRIM(name)//"_"//TRIM(int2str(ind_b))
272        ind_b=ind_b_in
273        ind_e=ind_b_in
274        single=.TRUE.
275     ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
276        name=TRIM(name)//"_"//TRIM(int2str(ind_e))
277        ind_b=ind_e_in
278        ind_e=ind_e_in
279        single=.TRUE.
280     ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
281        ind_b=ind_b_in
282        ind_e=ind_e_in
283        single=.FALSE.
284     ELSE
285        ind_b=1
286        ind_e=ndomain
287        single=.FALSE.
288     ENDIF
289     
290     NbField=NbField+1
291     FieldName(NbField)=TRIM(ADJUSTL(name))
292     FieldIndex(NbField)=1
293     
294!     PRINT *, 'Creating header for writefield : ', TRIM(FieldName(NbField)), ndomain, ind_b, ind_e ! FIXME
295!     PRINT *, mesh_glo%ndomain, SIZE(mesh_glo%primal_own)
296     SELECT CASE(Field(ind_b)%field_type)
297     CASE(field_T)
298!        PRINT *, '    field_type == field_T' ! FIXME
299        IF(single) THEN ! include halos
300           cells => mesh_glo%primal_all
301        ELSE
302           cells => mesh_glo%primal_own
303        END IF
304     CASE(field_Z)
305!        PRINT *, '    field_type == field_Z' ! FIXME
306        IF(single) THEN ! include halos
307           cells => mesh_glo%dual_all
308        ELSE
309           cells => mesh_glo%dual_own
310        END IF
311     END SELECT
312     
313!     PRINT *, 'writefield : ind_b, ind_e :', ind_b, ind_e
314
315     ncell=0
316     DO ind=ind_b,ind_e
317        nvert=SIZE(cells(ind)%bnds_lon,1)
318        ncell = ncell + cells(ind)%ncell
319     END DO
320!     PRINT *, 'writefield : found ',ncell, ' cells.' ! FIXME
321     
322     status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
323     FieldId(NbField)=ncid
324     status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
325     
326     level_size=0
327     SELECT CASE(Field(ind_b)%ndim)
328     CASE(2)
329        FieldVarId(NbField)%size=1
330        ALLOCATE(FieldVarId(NbField)%nc_id(1))
331     CASE(3)
332        FieldVarId(NbField)%size=1
333        ALLOCATE(FieldVarId(NbField)%nc_id(1))
334        status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
335        level_size=size(field(ind_b)%rval3d,2)
336     CASE(4)
337        FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
338        ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
339        status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
340        level_size=size(field(ind_b)%rval4d,2)
341        !          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
342     END SELECT
343     
344     PRINT*,"create_header_gen : LEVEL_SIZE=",level_size
345     IF (level_size>0) THEN
346        status = NF90_DEF_VAR(ncid,'lev',NF90_DOUBLE,(/ dim3id /),levId)
347        status = NF90_PUT_ATT(ncid,levId,"axis","Z")
348     ENDIF
349     
350     SELECT CASE(Field(ind_b)%field_type)
351     CASE(field_T)
352        status = NF90_DEF_DIM(ncid,'cell_i',ncell,ncellId)
353        status = NF90_DEF_DIM(ncid,'nvert_i',nvert,nvertid)
354        status = NF90_DEF_VAR(ncid,'lon_i',NF90_DOUBLE,(/ ncellId /),lonId)
355        status = NF90_DEF_VAR(ncid,'lat_i',NF90_DOUBLE,(/ ncellId /),latId)
356        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon_i")
357        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat_i")
358        status = NF90_DEF_VAR(ncid,'bounds_lon_i',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
359        status = NF90_DEF_VAR(ncid,'bounds_lat_i',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
360     CASE(field_Z)
361        status = NF90_DEF_DIM(ncid,'cell_v',ncell,ncellId)
362        status = NF90_DEF_DIM(ncid,'nvert_v',nvert,nvertid)
363        status = NF90_DEF_VAR(ncid,'lon_v',NF90_DOUBLE,(/ ncellId /),lonId)
364        status = NF90_DEF_VAR(ncid,'lat_v',NF90_DOUBLE,(/ ncellId /),latId)
365        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon_v")
366        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat_v")
367        status = NF90_DEF_VAR(ncid,'bounds_lon_v',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
368        status = NF90_DEF_VAR(ncid,'bounds_lat_v',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
369     END SELECT
370     
371     status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
372     status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
373     status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
374     status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
375     
376     SELECT CASE(Field(ind_b)%ndim)
377     CASE(2)
378        IF (once) THEN
379           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, &
380                (/ ncellId /),FieldVarId(NbField)%nc_id(1))
381        ELSE
382           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, &
383                (/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
384        END IF
385        CALL put_att_coordinates(1)
386     CASE(3)
387        IF (once) THEN
388           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, &
389                (/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(1))
390        ELSE
391           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, &
392                (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
393        END IF
394        CALL put_att_coordinates(1)
395     CASE(4)
396        DO i=1,FieldVarId(NbField)%size
397           IF (once) THEN
398              status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id /),  &
399                   FieldVarId(NbField)%nc_id(i))
400           ELSE
401              status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /),  &
402                   FieldVarId(NbField)%nc_id(i))
403           END IF
404           CALL put_att_coordinates(i)
405           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon_i lat_i")
406        END DO
407     END SELECT
408     
409     status = NF90_ENDDEF(ncid)     
410     
411     if (level_size>0) status = NF90_PUT_VAR(ncid,levId,(/ (l,l=1,level_size) /))
412     
413     ALLOCATE(lon(ncell),lat(ncell))
414     ALLOCATE(bounds_lon(nvert,ncell), bounds_lat(nvert,ncell))
415     n_begin=0 
416     DO ind=ind_b,ind_e
417        n_end = n_begin + cells(ind)%ncell
418!        PRINT *, 'create_header_gen ', n_begin, n_end, SHAPE(cells(ind)%bnds_lon), SHAPE(bounds_lon)
419        lat(n_begin+1:n_end) = cells(ind)%lat(:) *180./Pi
420        lon(n_begin+1:n_end) = cells(ind)%lon(:) *180./Pi
421        bounds_lon(:,n_begin+1:n_end) = cells(ind)%bnds_lon(:,:) *180./Pi
422        bounds_lat(:,n_begin+1:n_end) = cells(ind)%bnds_lat(:,:) *180./Pi
423        n_begin = n_end
424     END DO
425     status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /))
426     status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /))
427     status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
428     status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
429     
430     DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
431     
432   CONTAINS
433     
434     SUBROUTINE put_att_coordinates(ind)
435       INTEGER :: ind
436       SELECT CASE(Field(ind_b)%field_type)
437       CASE(field_T)
438          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(ind), &
439               "coordinates","lon_i lat_i")
440       CASE(field_Z)
441          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(ind), &
442               "coordinates","lon_v lat_v")     
443       END SELECT
444     END SUBROUTINE put_att_coordinates
445     
446   END SUBROUTINE Create_Header_gen
447   
448   SUBROUTINE Close_files
449   USE netcdf_mod
450     INTEGER :: i,k,status
451!$OMP MASTER     
452     DO i=1,NbField
453         status=NF90_CLOSE(FieldId(i))
454     ENDDO
455!$OMP END MASTER
456   END SUBROUTINE  Close_files
457     
458     
459end module write_field_mod
Note: See TracBrowser for help on using the repository browser.