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

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

Merge advection scheme from sarvesh in standard version

YM
File size: 20.0 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,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /))
157          DEALLOCATE(field_val2d)
158        ELSE IF (field(ind)%ndim==3) THEN
159          ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
160          n=0
161          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
162            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
163              k=d%iim*(j-1)+i
164              IF (d%own(i,j) .OR. single) THEN
165                n=n+1
166                Field_val3d(n,:)=field(ind)%rval3d(k,:)
167              ENDIF
168            ENDDO
169          ENDDO
170           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   &
171                               count=(/n,size(field(1)%rval3d,2),1 /))
172          DEALLOCATE(field_val3d)
173        ELSE IF (field(1)%ndim==4) THEN
174
175          DO q=1,FieldVarId(index)%size
176         
177            ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
178            n=0
179            DO j=d%jj_begin-halo_size,d%jj_end+halo_size
180              DO i=d%ii_begin-halo_size,d%ii_end+halo_size
181                k=d%iim*(j-1)+i
182                IF (d%own(i,j) .OR. single) THEN
183                  n=n+1
184                  Field_val3d(n,:)=field(ind)%rval4d(k,:,q)
185                ENDIF
186              ENDDO
187            ENDDO
188            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   &
189                               count=(/n,size(field(1)%rval4d,2),1 /))
190            DEALLOCATE(field_val3d)
191          ENDDO         
192        ENDIF
193       
194!        PRINT *,NF90_STRERROR(status)
195        ncell=ncell+n
196
197     ENDDO
198     
199     ELSE IF (Field(ind_b)%field_type==field_Z) THEN
200        ncell=1
201        n=0
202        DO ind=ind_b,ind_e
203          d=>domain(ind)
204          CALL swap_geometry(ind)
205          CALL swap_dimensions(ind)
206 
207          n=0
208          DO j=jj_begin+1,jj_end
209            DO i=ii_begin,ii_end-1
210              n=n+1
211            ENDDO
212          ENDDO
213
214          DO j=jj_begin,jj_end-1
215            DO i=ii_begin+1,ii_end
216              n=n+1
217            ENDDO
218          ENDDO
219
220        IF (field(ind)%ndim==2) THEN
221          ALLOCATE(Field_val2d(n))
222
223          n=0
224          DO j=jj_begin+1,jj_end
225            DO i=ii_begin,ii_end-1
226              n=n+1
227              k=iim*(j-1)+i
228              Field_val2d(n)=field(ind)%rval2d(k+z_down)
229            ENDDO
230          ENDDO
231
232          DO j=jj_begin,jj_end-1
233            DO i=ii_begin+1,ii_end
234              n=n+1
235              k=iim*(j-1)+i
236              Field_val2d(n)=field(ind)%rval2d(k+z_up)
237            ENDDO
238          ENDDO
239
240          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /))
241          DEALLOCATE(field_val2d)
242
243        ELSE IF (field(ind)%ndim==3) THEN
244          ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
245          n=0
246          DO j=jj_begin+1,jj_end
247            DO i=ii_begin,ii_end-1
248              n=n+1
249              k=iim*(j-1)+i
250              Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:)
251            ENDDO
252          ENDDO
253
254          DO j=jj_begin,jj_end-1
255            DO i=ii_begin+1,ii_end
256              n=n+1
257              k=iim*(j-1)+i
258              Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:)
259            ENDDO
260          ENDDO
261           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   &
262                               count=(/n,size(field(1)%rval3d,2),1 /))
263          DEALLOCATE(field_val3d)
264        ELSE IF (field(1)%ndim==4) THEN
265
266          DO q=1,FieldVarId(index)%size
267            ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
268            n=0
269            DO j=jj_begin+1,jj_end
270              DO i=ii_begin,ii_end-1
271                n=n+1
272                k=iim*(j-1)+i
273                Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q)
274              ENDDO
275            ENDDO
276
277            DO j=jj_begin,jj_end-1
278              DO i=ii_begin+1,ii_end
279                n=n+1
280                k=iim*(j-1)+i
281                Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q)
282              ENDDO
283            ENDDO
284
285            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,1,FieldIndex(Index) /),   &
286                                count=(/n,size(field(1)%rval4d,2),1 /))
287            DEALLOCATE(field_val3d)
288          ENDDO
289        ENDIF
290       
291!        PRINT *,NF90_STRERROR(status)
292        ncell=ncell+n
293
294     ENDDO
295     
296     ENDIF
297     status=NF90_SYNC(FieldId(Index))
298     
299    END SUBROUTINE Writefield
300     
301           
302    SUBROUTINE Create_header(name,field,nind)
303    USE netcdf
304    USE field_mod
305    USE domain_mod
306    USE spherical_geom_mod
307    USE dimensions
308    USE geometry
309    IMPLICIT NONE
310      CHARACTER(LEN=*) :: name
311      TYPE(t_field),POINTER :: field(:)
312      INTEGER,OPTIONAL,INTENT(IN) :: nind
313      INTEGER :: ncell
314      INTEGER :: nvert
315      REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
316      TYPE(t_domain),POINTER :: d
317      INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
318      INTEGER :: dim3id,dim4id
319      INTEGER :: status
320      INTEGER :: ind,i,j,k,n,q
321      INTEGER :: iie,jje,iin,jjn
322      INTEGER :: ind_b,ind_e
323      INTEGER :: halo_size
324      LOGICAL :: single 
325      INTEGER :: nij
326         
327      NbField=NbField+1
328      FieldName(NbField)=TRIM(ADJUSTL(name))
329      FieldIndex(NbField)=1
330     
331      IF (PRESENT(nind)) THEN
332        ind_b=nind
333        ind_e=nind
334        halo_size=1
335        single=.TRUE.
336      ELSE
337        ind_b=1
338        ind_e=ndomain
339        halo_size=0
340        single=.FALSE.
341      ENDIF
342     
343      ncell=0
344     
345      IF (Field(ind_b)%field_type==field_T) THEN
346        nvert=6
347       
348        DO ind=ind_b,ind_e
349          d=>domain(ind)
350       
351          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
352            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
353              IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1
354            ENDDO
355          ENDDO
356
357        END DO
358     
359        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
360        FieldId(NbField)=ncid
361        status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId)
362        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
363
364        IF (Field(ind_b)%ndim==2)  THEN
365          FieldVarId(NbField)%size=1
366          ALLOCATE(FieldVarId(NbField)%nc_id(1))
367        ELSE IF (Field(ind_b)%ndim==3)  THEN
368          FieldVarId(NbField)%size=1
369          ALLOCATE(FieldVarId(NbField)%nc_id(1))
370          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id)
371        ELSE IF (Field(1)%ndim==4) THEN
372          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
373          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
374          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id)
375!          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
376        ENDIF
377     
378        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
379     
380        status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
381        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
382        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
383        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
384        status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
385        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
386        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
387        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
388        status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
389        status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
390
391        IF (Field(ind_b)%ndim==2) THEN
392          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
393          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
394        ELSE IF (Field(ind_b)%ndim==3) THEN
395          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
396          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
397        ELSE IF (Field(ind_b)%ndim==4) THEN
398          DO i=1,FieldVarId(NbField)%size
399            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(i))
400            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat")
401          ENDDO       
402        ENDIF
403 
404     
405        status = NF90_ENDDEF(ncid)     
406
407        ncell=1
408        DO ind=ind_b,ind_e
409          d=>domain(ind)
410 
411          n=0
412          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
413            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
414              IF (single .OR. d%own(i,j)) n=n+1
415            ENDDO
416          ENDDO
417         
418         ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
419         
420          n=0 
421          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
422            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
423                IF (d%own(i,j) .OR. single) THEN
424                n=n+1
425                CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n))
426                lon(n)=lon(n)*180/Pi
427                lat(n)=lat(n)*180/Pi
428                DO k=0,5
429                  CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n))
430                  bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
431                  bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
432                ENDDO
433              ENDIF
434            ENDDO
435          ENDDO
436          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /))
437          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /))
438          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
439          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
440 
441          ncell=ncell+n
442          DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
443      END DO
444
445    ELSE IF (Field(ind_b)%field_type==field_Z) THEN
446        nvert=3
447        DO ind=ind_b,ind_e
448          d=>domain(ind)
449       
450          DO j=d%jj_begin+1,d%jj_end
451            DO i=d%ii_begin,d%ii_end-1
452              ncell=ncell+1
453            ENDDO
454          ENDDO
455
456          DO j=d%jj_begin,d%jj_end-1
457            DO i=d%ii_begin+1,d%ii_end
458              ncell=ncell+1
459            ENDDO
460          ENDDO
461
462        END DO
463     
464        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
465        FieldId(NbField)=ncid
466        status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId)
467        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
468
469        IF (Field(ind_b)%ndim==2)  THEN
470          FieldVarId(NbField)%size=1
471          ALLOCATE(FieldVarId(NbField)%nc_id(1))
472        ELSE IF (Field(ind_b)%ndim==3)  THEN
473          FieldVarId(NbField)%size=1
474          ALLOCATE(FieldVarId(NbField)%nc_id(1))
475          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id)
476        ELSE IF (Field(1)%ndim==4) THEN
477          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
478          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
479          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id)
480        ENDIF
481
482
483     
484        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
485     
486        status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
487        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
488        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
489        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
490        status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
491        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
492        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
493        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
494        status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
495        status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
496
497
498        IF (Field(ind_b)%ndim==2) THEN
499          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
500          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
501        ELSE IF (Field(ind_b)%ndim==3) THEN
502          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,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==4) THEN
505          DO q=1,FieldVarId(NbField)%size
506            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q))
507            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat")
508          ENDDO       
509        ENDIF
510       
511        status = NF90_ENDDEF(ncid)     
512
513        ncell=1
514        DO ind=ind_b,ind_e
515          d=>domain(ind)
516          CALL swap_geometry(ind)
517          CALL swap_dimensions(ind)
518 
519          n=0
520          DO j=jj_begin+1,jj_end
521            DO i=ii_begin,ii_end-1
522              n=n+1
523            ENDDO
524          ENDDO
525
526          DO j=jj_begin,jj_end-1
527            DO i=ii_begin+1,ii_end
528              n=n+1
529            ENDDO
530          ENDDO
531
532         ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
533         
534          n=0 
535       
536          DO j=jj_begin+1,jj_end
537            DO i=ii_begin,ii_end-1
538              nij=(j-1)*iim+i
539              n=n+1
540              CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n))
541              lon(n)=lon(n)*180/Pi
542 !             IF (lon(n)<0) lon(n)=lon(n)+360
543              lat(n)=lat(n)*180/Pi
544              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
545              CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
546              CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
547
548              DO k=0,2
549                bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
550                bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
551!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
552              ENDDO
553            ENDDO
554          ENDDO
555
556          DO j=jj_begin,jj_end-1
557            DO i=ii_begin+1,ii_end
558              nij=(j-1)*iim+i
559              n=n+1
560              CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
561              lon(n)=lon(n)*180/Pi
562!              IF (lon(n)<0) lon(n)=lon(n)+360
563              lat(n)=lat(n)*180/Pi
564              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
565              CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
566              CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
567
568              DO k=0,2
569                bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
570                bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
571!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
572              ENDDO
573            ENDDO
574          ENDDO
575         
576         
577          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /))
578          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /))
579          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
580          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
581 
582          ncell=ncell+n
583          DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
584      END DO         
585    ENDIF
586
587
588     
589!      status = NF90_CLOSE(ncid)
590
591    end subroutine Create_Header
592   
593 
594  function int2str(int)
595    implicit none
596    integer, parameter :: MaxLen=10
597    integer,intent(in) :: int
598    character(len=MaxLen) :: int2str
599    logical :: flag
600    integer :: i
601    flag=.true.
602   
603    i=int
604   
605    int2str=''
606    do while (flag)
607      int2str=CHAR(MOD(i,10)+48)//int2str
608      i=i/10
609      if (i==0) flag=.false.
610    enddo
611  end function int2str
612
613end module write_field
614 
Note: See TracBrowser for help on using the repository browser.