module write_field USE genmod implicit none integer, parameter :: MaxWriteField = 1000 integer, dimension(MaxWriteField),save :: FieldId integer, dimension(MaxWriteField),save :: FieldVarId integer, dimension(MaxWriteField),save :: FieldIndex character(len=255), dimension(MaxWriteField) :: FieldName integer,save :: NbField = 0 contains function GetFieldIndex(name) implicit none integer :: GetFieldindex character(len=*) :: name character(len=255) :: TrueName integer :: i TrueName=TRIM(ADJUSTL(name)) GetFieldIndex=-1 do i=1,NbField if (TrueName==FieldName(i)) then GetFieldIndex=i exit endif enddo end function GetFieldIndex subroutine WriteField_gen(name,Field,dimx,dimy,dimz) implicit none ! include 'netcdf.inc' character(len=*) :: name integer :: dimx,dimy,dimz real,dimension(dimx,dimy,dimz) :: Field integer,dimension(dimx*dimy*dimz) :: ndex integer :: status integer :: index integer :: start(4) integer :: count(4) Index=GetFieldIndex(name) if (Index==-1) then ! call CreateNewField(name,dimx,dimy,dimz) Index=GetFieldIndex(name) else FieldIndex(Index)=FieldIndex(Index)+1. endif start(1)=1 start(2)=1 start(3)=1 start(4)=FieldIndex(Index) count(1)=dimx count(2)=dimy count(3)=dimz count(4)=1 ! status = NF_PUT_VARA_DOUBLE(FieldId(Index),FieldVarId(Index),start,count,Field) ! status = NF_SYNC(FieldId(Index)) end subroutine WriteField_gen SUBROUTINE Writefield(name_in,field,nind) USE netcdf USE domain_mod use field_mod USE dimensions USE geometry IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: name_in TYPE(t_field),POINTER :: field(:) INTEGER,OPTIONAL,INTENT(IN) :: nind REAL(r8),ALLOCATABLE :: field_val2d(:) REAL(r8),ALLOCATABLE :: field_val3d(:,:) REAL(r8),ALLOCATABLE :: field_val4d(:,:,:) TYPE(t_domain),POINTER :: d INTEGER :: Index INTEGER :: ind,i,j,k,n,ncell INTEGER :: iie,jje,iin,jjn INTEGER :: status CHARACTER(len=255) :: name CHARACTER(len=255) :: str_ind INTEGER :: ind_b,ind_e INTEGER :: halo_size LOGICAL :: single name=TRIM(ADJUSTL(name_in)) IF (PRESENT(nind)) THEN name=TRIM(name)//"_"//TRIM(int2str(nind)) PRINT *,"NAME",nind,int2str(nind),name ind_b=nind ind_e=nind halo_size=1 single=.TRUE. ELSE ind_b=1 ind_e=ndomain halo_size=0 single=.FALSE. ENDIF Index=GetFieldIndex(name) if (Index==-1) then call create_header(name,field,nind) Index=GetFieldIndex(name) else FieldIndex(Index)=FieldIndex(Index)+1. endif IF (Field(ind_b)%field_type==field_T) THEN ncell=1 DO ind=ind_b,ind_e d=>domain(ind) IF (Field(ind)%field_type/=field_T) THEN PRINT *,"Writefield, grille non geree" RETURN ENDIF n=0 DO j=d%jj_begin-halo_size,d%jj_end+halo_size DO i=d%ii_begin-halo_size,d%ii_end+halo_size IF (d%own(i,j) .OR. single) n=n+1 ENDDO ENDDO IF (field(ind)%ndim==2) THEN ALLOCATE(Field_val2d(n)) n=0 DO j=d%jj_begin-halo_size,d%jj_end+halo_size DO i=d%ii_begin-halo_size,d%ii_end+halo_size k=d%iim*(j-1)+i IF (d%own(i,j) .OR. single) THEN n=n+1 Field_val2d(n)=field(ind)%rval2d(k) ENDIF ENDDO ENDDO status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index),Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) DEALLOCATE(field_val2d) ELSE IF (field(ind)%ndim==3) THEN ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2))) n=0 DO j=d%jj_begin-halo_size,d%jj_end+halo_size DO i=d%ii_begin-halo_size,d%ii_end+halo_size k=d%iim*(j-1)+i IF (d%own(i,j) .OR. single) THEN n=n+1 Field_val3d(n,:)=field(ind)%rval3d(k,:) ENDIF ENDDO ENDDO status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /), & count=(/n,size(field(1)%rval3d,2),1 /)) DEALLOCATE(field_val3d) ELSE IF (field(1)%ndim==4) THEN ALLOCATE(Field_val4d(n,size(field(ind)%rval4d,2),size(field(ind)%rval4d,3))) n=0 DO j=d%jj_begin-halo_size,d%jj_end+halo_size DO i=d%ii_begin-halo_size,d%ii_end+halo_size k=d%iim*(j-1)+i IF (d%own(i,j) .OR. single) THEN n=n+1 Field_val4d(n,:,:)=field(ind)%rval4d(k,:,:) ENDIF ENDDO ENDDO status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index),Field_val4d,start=(/ ncell,1,1,FieldIndex(Index) /), & count=(/n,size(field(1)%rval4d,2),size(field(1)%rval4d,3),1 /)) DEALLOCATE(field_val4d) ENDIF ! PRINT *,NF90_STRERROR(status) ncell=ncell+n ENDDO ELSE IF (Field(ind_b)%field_type==field_Z) THEN ncell=1 n=0 DO ind=ind_b,ind_e d=>domain(ind) CALL swap_geometry(ind) CALL swap_dimensions(ind) n=0 DO j=jj_begin+1,jj_end DO i=ii_begin,ii_end-1 n=n+1 ENDDO ENDDO DO j=jj_begin,jj_end-1 DO i=ii_begin+1,ii_end n=n+1 ENDDO ENDDO IF (field(ind)%ndim==2) THEN ALLOCATE(Field_val2d(n)) n=0 DO j=jj_begin+1,jj_end DO i=ii_begin,ii_end-1 n=n+1 k=iim*(j-1)+i Field_val2d(n)=field(ind)%rval2d(k+z_down) ENDDO ENDDO DO j=jj_begin,jj_end-1 DO i=ii_begin+1,ii_end n=n+1 k=iim*(j-1)+i Field_val2d(n)=field(ind)%rval2d(k+z_up) ENDDO ENDDO status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index),Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) DEALLOCATE(field_val2d) ELSE IF (field(ind)%ndim==3) THEN ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2))) n=0 DO j=jj_begin+1,jj_end DO i=ii_begin,ii_end-1 n=n+1 k=iim*(j-1)+i Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:) ENDDO ENDDO DO j=jj_begin,jj_end-1 DO i=ii_begin+1,ii_end n=n+1 k=iim*(j-1)+i Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:) ENDDO ENDDO status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /), & count=(/n,size(field(1)%rval3d,2),1 /)) DEALLOCATE(field_val3d) ELSE IF (field(1)%ndim==4) THEN ALLOCATE(Field_val4d(n,size(field(ind)%rval4d,2),size(field(ind)%rval4d,3))) n=0 DO j=jj_begin+1,jj_end DO i=ii_begin,ii_end-1 n=n+1 k=iim*(j-1)+i Field_val4d(n,:,:)=field(ind)%rval4d(k+z_down,:,:) ENDDO ENDDO DO j=jj_begin,jj_end-1 DO i=ii_begin+1,ii_end n=n+1 k=iim*(j-1)+i Field_val4d(n,:,:)=field(ind)%rval4d(k+z_up,:,:) ENDDO ENDDO status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index),Field_val4d,start=(/ ncell,1,1,FieldIndex(Index) /), & count=(/n,size(field(1)%rval4d,2),size(field(1)%rval4d,3),1 /)) DEALLOCATE(field_val4d) ENDIF ! PRINT *,NF90_STRERROR(status) ncell=ncell+n ENDDO ENDIF status=NF90_SYNC(FieldId(Index)) END SUBROUTINE Writefield SUBROUTINE Create_header(name,field,nind) USE netcdf USE field_mod USE domain_mod USE spherical_geom_mod USE dimensions USE geometry IMPLICIT NONE CHARACTER(LEN=*) :: name TYPE(t_field),POINTER :: field(:) INTEGER,OPTIONAL,INTENT(IN) :: nind INTEGER :: ncell INTEGER :: nvert REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:) TYPE(t_domain),POINTER :: d INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId INTEGER :: dim3id,dim4id INTEGER :: status INTEGER :: ind,i,j,k,n INTEGER :: iie,jje,iin,jjn INTEGER :: ind_b,ind_e INTEGER :: halo_size LOGICAL :: single INTEGER :: nij NbField=NbField+1 FieldName(NbField)=TRIM(ADJUSTL(name)) FieldIndex(NbField)=1 IF (PRESENT(nind)) THEN ind_b=nind ind_e=nind halo_size=1 single=.TRUE. ELSE ind_b=1 ind_e=ndomain halo_size=0 single=.FALSE. ENDIF ncell=0 IF (Field(ind_b)%field_type==field_T) THEN nvert=6 DO ind=ind_b,ind_e d=>domain(ind) DO j=d%jj_begin-halo_size,d%jj_end+halo_size DO i=d%ii_begin-halo_size,d%ii_end+halo_size IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1 ENDDO ENDDO END DO status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) FieldId(NbField)=ncid status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) IF (Field(ind_b)%ndim==3) THEN status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) ELSE IF (Field(1)%ndim==4) THEN status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id) ENDIF status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) IF (Field(ind_b)%ndim==2) THEN status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)) ELSE IF (Field(ind_b)%ndim==3) THEN status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)) ELSE IF (Field(ind_b)%ndim==4) THEN status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,dim4id,timeId /),FieldVarId(NbField)) ENDIF status = NF90_PUT_ATT(ncid,FieldVarId(NbField),"coordinates","lon lat") status = NF90_ENDDEF(ncid) ncell=1 DO ind=ind_b,ind_e d=>domain(ind) n=0 DO j=d%jj_begin-halo_size,d%jj_end+halo_size DO i=d%ii_begin-halo_size,d%ii_end+halo_size IF (single .OR. d%own(i,j)) n=n+1 ENDDO ENDDO ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n)) n=0 DO j=d%jj_begin-halo_size,d%jj_end+halo_size DO i=d%ii_begin-halo_size,d%ii_end+halo_size IF (d%own(i,j) .OR. single) THEN n=n+1 CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n)) lon(n)=lon(n)*180/Pi lat(n)=lat(n)*180/Pi DO k=0,5 CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n)) bounds_lat(k,n)=bounds_lat(k,n)*180/Pi bounds_lon(k,n)=bounds_lon(k,n)*180/Pi ENDDO ENDIF ENDDO ENDDO status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /)) status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /)) status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) ncell=ncell+n DEALLOCATE(lon,lat,bounds_lon,bounds_lat) END DO ELSE IF (Field(ind_b)%field_type==field_Z) THEN nvert=3 DO ind=ind_b,ind_e d=>domain(ind) DO j=d%jj_begin+1,d%jj_end DO i=d%ii_begin,d%ii_end-1 ncell=ncell+1 ENDDO ENDDO DO j=d%jj_begin,d%jj_end-1 DO i=d%ii_begin+1,d%ii_end ncell=ncell+1 ENDDO ENDDO END DO status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) FieldId(NbField)=ncid status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) IF (Field(ind_b)%ndim==3) THEN status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) ELSE IF (Field(1)%ndim==4) THEN status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id) ENDIF status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) IF (Field(ind_b)%ndim==2) THEN status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)) ELSE IF (Field(ind_b)%ndim==3) THEN status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)) ELSE IF (Field(ind_b)%ndim==4) THEN status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,dim4id,timeId /),FieldVarId(NbField)) ENDIF status = NF90_PUT_ATT(ncid,FieldVarId(NbField),"coordinates","lon lat") status = NF90_ENDDEF(ncid) ncell=1 DO ind=ind_b,ind_e d=>domain(ind) CALL swap_geometry(ind) CALL swap_dimensions(ind) n=0 DO j=jj_begin+1,jj_end DO i=ii_begin,ii_end-1 n=n+1 ENDDO ENDDO DO j=jj_begin,jj_end-1 DO i=ii_begin+1,ii_end n=n+1 ENDDO ENDDO ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n)) n=0 DO j=jj_begin+1,jj_end DO i=ii_begin,ii_end-1 nij=(j-1)*iim+i n=n+1 CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n)) lon(n)=lon(n)*180/Pi ! IF (lon(n)<0) lon(n)=lon(n)+360 lat(n)=lat(n)*180/Pi CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) DO k=0,2 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi bounds_lon(k,n)=bounds_lon(k,n)*180/Pi ! IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 ENDDO ENDDO ENDDO DO j=jj_begin,jj_end-1 DO i=ii_begin+1,ii_end nij=(j-1)*iim+i n=n+1 CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n)) lon(n)=lon(n)*180/Pi ! IF (lon(n)<0) lon(n)=lon(n)+360 lat(n)=lat(n)*180/Pi CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) DO k=0,2 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi bounds_lon(k,n)=bounds_lon(k,n)*180/Pi ! IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 ENDDO ENDDO ENDDO status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /)) status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /)) status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) ncell=ncell+n DEALLOCATE(lon,lat,bounds_lon,bounds_lat) END DO ENDIF ! status = NF90_CLOSE(ncid) end subroutine Create_Header function int2str(int) implicit none integer, parameter :: MaxLen=10 integer,intent(in) :: int character(len=MaxLen) :: int2str logical :: flag integer :: i flag=.true. i=int int2str='' do while (flag) int2str=CHAR(MOD(i,10)+48)//int2str i=i/10 if (i==0) flag=.false. enddo end function int2str end module write_field