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

Last change on this file since 21 was 21, checked in by ymipsl, 12 years ago

correction for compiling with gfortran (line too long)
improvement for splitting domain
Call twice transfert request for field u is no longer necessary

YM
File size: 20.2 KB
Line 
1module write_field
2USE genmod
3implicit none
4
5  TYPE ncvar
6    INTEGER :: size
7    INTEGER,POINTER :: nc_id(:)
8  END TYPE ncvar
9
10  integer, parameter :: MaxWriteField = 1000
11  integer, dimension(MaxWriteField),save :: FieldId
12  TYPE(ncvar), dimension(MaxWriteField),save :: FieldVarId
13  integer, dimension(MaxWriteField),save :: FieldIndex
14  character(len=255), dimension(MaxWriteField) ::  FieldName 
15   
16  integer,save :: NbField = 0
17 
18  contains
19 
20    function GetFieldIndex(name)
21    implicit none
22      integer          :: GetFieldindex
23      character(len=*) :: name
24   
25      character(len=255) :: TrueName
26      integer            :: i
27       
28     
29      TrueName=TRIM(ADJUSTL(name))
30   
31      GetFieldIndex=-1
32      do i=1,NbField
33        if (TrueName==FieldName(i)) then
34          GetFieldIndex=i
35          exit
36        endif
37      enddo
38    end function GetFieldIndex
39 
40
41     
42    subroutine WriteField_gen(name,Field,dimx,dimy,dimz)
43    implicit none
44!    include 'netcdf.inc'
45      character(len=*) :: name
46      integer :: dimx,dimy,dimz
47      real,dimension(dimx,dimy,dimz) :: Field
48      integer,dimension(dimx*dimy*dimz) :: ndex
49      integer :: status
50      integer :: index
51      integer :: start(4)
52      integer :: count(4)
53     
54           
55      Index=GetFieldIndex(name)
56      if (Index==-1) then
57!        call CreateNewField(name,dimx,dimy,dimz)
58        Index=GetFieldIndex(name)
59      else
60        FieldIndex(Index)=FieldIndex(Index)+1.
61      endif
62     
63      start(1)=1
64      start(2)=1
65      start(3)=1
66      start(4)=FieldIndex(Index)
67
68      count(1)=dimx
69      count(2)=dimy
70      count(3)=dimz
71      count(4)=1
72
73!      status = NF_PUT_VARA_DOUBLE(FieldId(Index),FieldVarId(Index),start,count,Field)
74!      status = NF_SYNC(FieldId(Index))
75     
76    end subroutine WriteField_gen
77
78
79    SUBROUTINE Writefield(name_in,field,nind)
80    USE netcdf
81    USE domain_mod
82    use field_mod
83    USE dimensions
84    USE geometry
85    IMPLICIT NONE
86      CHARACTER(LEN=*),INTENT(IN) :: name_in
87      TYPE(t_field),POINTER :: field(:)
88      INTEGER,OPTIONAL,INTENT(IN) :: nind
89      REAL(r8),ALLOCATABLE :: field_val2d(:)
90      REAL(r8),ALLOCATABLE :: field_val3d(:,:)
91      REAL(r8),ALLOCATABLE :: field_val4d(:,:,:)
92      TYPE(t_domain),POINTER :: d
93      INTEGER :: Index
94      INTEGER :: ind,i,j,k,n,ncell,q
95      INTEGER :: iie,jje,iin,jjn
96      INTEGER :: status
97      CHARACTER(len=255) :: name
98      CHARACTER(len=255) :: str_ind
99      INTEGER :: ind_b,ind_e
100      INTEGER :: halo_size
101      LOGICAL :: single
102     
103     
104      name=TRIM(ADJUSTL(name_in))
105
106      IF (PRESENT(nind)) THEN
107        name=TRIM(name)//"_"//TRIM(int2str(nind))
108        PRINT *,"NAME",nind,int2str(nind),name
109        ind_b=nind
110        ind_e=nind
111        halo_size=1
112        single=.TRUE.
113      ELSE
114        ind_b=1
115        ind_e=ndomain
116        halo_size=0
117        single=.FALSE.
118      ENDIF     
119
120      Index=GetFieldIndex(name)
121      if (Index==-1) then
122        call create_header(name,field,nind)
123        Index=GetFieldIndex(name)
124      else
125        FieldIndex(Index)=FieldIndex(Index)+1.
126      endif
127     
128      IF (Field(ind_b)%field_type==field_T) THEN
129        ncell=1
130        DO ind=ind_b,ind_e
131          d=>domain(ind)
132          IF (Field(ind)%field_type/=field_T) THEN
133            PRINT *,"Writefield, grille non geree"
134            RETURN
135          ENDIF
136
137        n=0
138        DO j=d%jj_begin-halo_size,d%jj_end+halo_size
139          DO i=d%ii_begin-halo_size,d%ii_end+halo_size
140            IF (d%own(i,j) .OR. single) n=n+1
141          ENDDO
142        ENDDO
143
144        IF (field(ind)%ndim==2) THEN
145          ALLOCATE(Field_val2d(n))
146          n=0
147          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
148            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
149              k=d%iim*(j-1)+i
150              IF (d%own(i,j) .OR. single) THEN
151                n=n+1
152                Field_val2d(n)=field(ind)%rval2d(k)
153              ENDIF
154            ENDDO
155          ENDDO
156          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  &
157                              start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /))
158          DEALLOCATE(field_val2d)
159        ELSE IF (field(ind)%ndim==3) THEN
160          ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
161          n=0
162          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
163            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
164              k=d%iim*(j-1)+i
165              IF (d%own(i,j) .OR. single) THEN
166                n=n+1
167                Field_val3d(n,:)=field(ind)%rval3d(k,:)
168              ENDIF
169            ENDDO
170          ENDDO
171           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   &
172                               count=(/n,size(field(1)%rval3d,2),1 /))
173          DEALLOCATE(field_val3d)
174        ELSE IF (field(1)%ndim==4) THEN
175
176          DO q=1,FieldVarId(index)%size
177         
178            ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
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_val3d(n,:)=field(ind)%rval4d(k,:,q)
186                ENDIF
187              ENDDO
188            ENDDO
189            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   &
190                               count=(/n,size(field(1)%rval4d,2),1 /))
191            DEALLOCATE(field_val3d)
192          ENDDO         
193        ENDIF
194       
195!        PRINT *,NF90_STRERROR(status)
196        ncell=ncell+n
197
198     ENDDO
199     
200     ELSE IF (Field(ind_b)%field_type==field_Z) THEN
201        ncell=1
202        n=0
203        DO ind=ind_b,ind_e
204          d=>domain(ind)
205          CALL swap_geometry(ind)
206          CALL swap_dimensions(ind)
207 
208          n=0
209          DO j=jj_begin+1,jj_end
210            DO i=ii_begin,ii_end-1
211              n=n+1
212            ENDDO
213          ENDDO
214
215          DO j=jj_begin,jj_end-1
216            DO i=ii_begin+1,ii_end
217              n=n+1
218            ENDDO
219          ENDDO
220
221        IF (field(ind)%ndim==2) THEN
222          ALLOCATE(Field_val2d(n))
223
224          n=0
225          DO j=jj_begin+1,jj_end
226            DO i=ii_begin,ii_end-1
227              n=n+1
228              k=iim*(j-1)+i
229              Field_val2d(n)=field(ind)%rval2d(k+z_down)
230            ENDDO
231          ENDDO
232
233          DO j=jj_begin,jj_end-1
234            DO i=ii_begin+1,ii_end
235              n=n+1
236              k=iim*(j-1)+i
237              Field_val2d(n)=field(ind)%rval2d(k+z_up)
238            ENDDO
239          ENDDO
240
241          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       &
242                              Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /))
243          DEALLOCATE(field_val2d)
244
245        ELSE IF (field(ind)%ndim==3) THEN
246          ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
247          n=0
248          DO j=jj_begin+1,jj_end
249            DO i=ii_begin,ii_end-1
250              n=n+1
251              k=iim*(j-1)+i
252              Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:)
253            ENDDO
254          ENDDO
255
256          DO j=jj_begin,jj_end-1
257            DO i=ii_begin+1,ii_end
258              n=n+1
259              k=iim*(j-1)+i
260              Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:)
261            ENDDO
262          ENDDO
263           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   &
264                               count=(/n,size(field(1)%rval3d,2),1 /))
265          DEALLOCATE(field_val3d)
266        ELSE IF (field(1)%ndim==4) THEN
267
268          DO q=1,FieldVarId(index)%size
269            ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
270            n=0
271            DO j=jj_begin+1,jj_end
272              DO i=ii_begin,ii_end-1
273                n=n+1
274                k=iim*(j-1)+i
275                Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q)
276              ENDDO
277            ENDDO
278
279            DO j=jj_begin,jj_end-1
280              DO i=ii_begin+1,ii_end
281                n=n+1
282                k=iim*(j-1)+i
283                Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q)
284              ENDDO
285            ENDDO
286
287            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,1,FieldIndex(Index) /),   &
288                                count=(/n,size(field(1)%rval4d,2),1 /))
289            DEALLOCATE(field_val3d)
290          ENDDO
291        ENDIF
292       
293!        PRINT *,NF90_STRERROR(status)
294        ncell=ncell+n
295
296     ENDDO
297     
298     ENDIF
299     status=NF90_SYNC(FieldId(Index))
300     
301    END SUBROUTINE Writefield
302     
303           
304    SUBROUTINE Create_header(name,field,nind)
305    USE netcdf
306    USE field_mod
307    USE domain_mod
308    USE spherical_geom_mod
309    USE dimensions
310    USE geometry
311    IMPLICIT NONE
312      CHARACTER(LEN=*) :: name
313      TYPE(t_field),POINTER :: field(:)
314      INTEGER,OPTIONAL,INTENT(IN) :: nind
315      INTEGER :: ncell
316      INTEGER :: nvert
317      REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
318      TYPE(t_domain),POINTER :: d
319      INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
320      INTEGER :: dim3id,dim4id
321      INTEGER :: status
322      INTEGER :: ind,i,j,k,n,q
323      INTEGER :: iie,jje,iin,jjn
324      INTEGER :: ind_b,ind_e
325      INTEGER :: halo_size
326      LOGICAL :: single 
327      INTEGER :: nij
328         
329      NbField=NbField+1
330      FieldName(NbField)=TRIM(ADJUSTL(name))
331      FieldIndex(NbField)=1
332     
333      IF (PRESENT(nind)) THEN
334        ind_b=nind
335        ind_e=nind
336        halo_size=1
337        single=.TRUE.
338      ELSE
339        ind_b=1
340        ind_e=ndomain
341        halo_size=0
342        single=.FALSE.
343      ENDIF
344     
345      ncell=0
346     
347      IF (Field(ind_b)%field_type==field_T) THEN
348        nvert=6
349       
350        DO ind=ind_b,ind_e
351          d=>domain(ind)
352       
353          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
354            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
355              IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1
356            ENDDO
357          ENDDO
358
359        END DO
360     
361        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
362        FieldId(NbField)=ncid
363        status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId)
364        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
365
366        IF (Field(ind_b)%ndim==2)  THEN
367          FieldVarId(NbField)%size=1
368          ALLOCATE(FieldVarId(NbField)%nc_id(1))
369        ELSE IF (Field(ind_b)%ndim==3)  THEN
370          FieldVarId(NbField)%size=1
371          ALLOCATE(FieldVarId(NbField)%nc_id(1))
372          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id)
373        ELSE IF (Field(1)%ndim==4) THEN
374          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
375          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
376          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id)
377!          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
378        ENDIF
379     
380        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
381     
382        status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
383        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
384        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
385        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
386        status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
387        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
388        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
389        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
390        status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
391        status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
392
393        IF (Field(ind_b)%ndim==2) THEN
394          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
395          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
396        ELSE IF (Field(ind_b)%ndim==3) THEN
397          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
398          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
399        ELSE IF (Field(ind_b)%ndim==4) THEN
400          DO i=1,FieldVarId(NbField)%size
401            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),  &
402                                  FieldVarId(NbField)%nc_id(i))
403            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat")
404          ENDDO       
405        ENDIF
406 
407     
408        status = NF90_ENDDEF(ncid)     
409
410        ncell=1
411        DO ind=ind_b,ind_e
412          d=>domain(ind)
413 
414          n=0
415          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
416            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
417              IF (single .OR. d%own(i,j)) n=n+1
418            ENDDO
419          ENDDO
420         
421         ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
422         
423          n=0 
424          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
425            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
426                IF (d%own(i,j) .OR. single) THEN
427                n=n+1
428                CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n))
429                lon(n)=lon(n)*180/Pi
430                lat(n)=lat(n)*180/Pi
431                DO k=0,5
432                  CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n))
433                  bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
434                  bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
435                ENDDO
436              ENDIF
437            ENDDO
438          ENDDO
439          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /))
440          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /))
441          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
442          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
443 
444          ncell=ncell+n
445          DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
446      END DO
447
448    ELSE IF (Field(ind_b)%field_type==field_Z) THEN
449        nvert=3
450        DO ind=ind_b,ind_e
451          d=>domain(ind)
452       
453          DO j=d%jj_begin+1,d%jj_end
454            DO i=d%ii_begin,d%ii_end-1
455              ncell=ncell+1
456            ENDDO
457          ENDDO
458
459          DO j=d%jj_begin,d%jj_end-1
460            DO i=d%ii_begin+1,d%ii_end
461              ncell=ncell+1
462            ENDDO
463          ENDDO
464
465        END DO
466     
467        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
468        FieldId(NbField)=ncid
469        status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId)
470        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
471
472        IF (Field(ind_b)%ndim==2)  THEN
473          FieldVarId(NbField)%size=1
474          ALLOCATE(FieldVarId(NbField)%nc_id(1))
475        ELSE IF (Field(ind_b)%ndim==3)  THEN
476          FieldVarId(NbField)%size=1
477          ALLOCATE(FieldVarId(NbField)%nc_id(1))
478          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id)
479        ELSE IF (Field(1)%ndim==4) THEN
480          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
481          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
482          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id)
483        ENDIF
484
485
486     
487        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
488     
489        status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
490        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
491        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
492        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
493        status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
494        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
495        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
496        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
497        status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
498        status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
499
500
501        IF (Field(ind_b)%ndim==2) THEN
502          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
503          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
504        ELSE IF (Field(ind_b)%ndim==3) THEN
505          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
506          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
507        ELSE IF (Field(ind_b)%ndim==4) THEN
508          DO q=1,FieldVarId(NbField)%size
509            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),NF90_DOUBLE,             &
510                                  (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q))
511            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat")
512          ENDDO       
513        ENDIF
514       
515        status = NF90_ENDDEF(ncid)     
516
517        ncell=1
518        DO ind=ind_b,ind_e
519          d=>domain(ind)
520          CALL swap_geometry(ind)
521          CALL swap_dimensions(ind)
522 
523          n=0
524          DO j=jj_begin+1,jj_end
525            DO i=ii_begin,ii_end-1
526              n=n+1
527            ENDDO
528          ENDDO
529
530          DO j=jj_begin,jj_end-1
531            DO i=ii_begin+1,ii_end
532              n=n+1
533            ENDDO
534          ENDDO
535
536         ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
537         
538          n=0 
539       
540          DO j=jj_begin+1,jj_end
541            DO i=ii_begin,ii_end-1
542              nij=(j-1)*iim+i
543              n=n+1
544              CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n))
545              lon(n)=lon(n)*180/Pi
546 !             IF (lon(n)<0) lon(n)=lon(n)+360
547              lat(n)=lat(n)*180/Pi
548              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
549              CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
550              CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
551
552              DO k=0,2
553                bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
554                bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
555!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
556              ENDDO
557            ENDDO
558          ENDDO
559
560          DO j=jj_begin,jj_end-1
561            DO i=ii_begin+1,ii_end
562              nij=(j-1)*iim+i
563              n=n+1
564              CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
565              lon(n)=lon(n)*180/Pi
566!              IF (lon(n)<0) lon(n)=lon(n)+360
567              lat(n)=lat(n)*180/Pi
568              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
569              CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
570              CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
571
572              DO k=0,2
573                bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
574                bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
575!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
576              ENDDO
577            ENDDO
578          ENDDO
579         
580         
581          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /))
582          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /))
583          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
584          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
585 
586          ncell=ncell+n
587          DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
588      END DO         
589    ENDIF
590
591
592     
593!      status = NF90_CLOSE(ncid)
594
595    end subroutine Create_Header
596   
597 
598  function int2str(int)
599    implicit none
600    integer, parameter :: MaxLen=10
601    integer,intent(in) :: int
602    character(len=MaxLen) :: int2str
603    logical :: flag
604    integer :: i
605    flag=.true.
606   
607    i=int
608   
609    int2str=''
610    do while (flag)
611      int2str=CHAR(MOD(i,10)+48)//int2str
612      i=i/10
613      if (i==0) flag=.false.
614    enddo
615  end function int2str
616
617end module write_field
618 
Note: See TracBrowser for help on using the repository browser.