module write_field_mpi_mod USE genmod USE write_field_vars_mod IMPLICIT NONE PRIVATE ! TYPE ncvar ! INTEGER :: size ! INTEGER,POINTER :: nc_id(:) ! INTEGER :: displ ! END TYPE ncvar ! ! INTEGER, PARAMETER :: MaxWriteField = 1000 ! INTEGER, DIMENSION(MaxWriteField),SAVE :: FieldId ! TYPE(ncvar), dimension(MaxWriteField),SAVE :: FieldVarId ! INTEGER, DIMENSION(MaxWriteField),SAVE :: FieldIndex ! CHARACTER(len=255), DIMENSION(MaxWriteField) :: FieldName ! ! INTEGER,SAVE :: NbField = 0 ! ! PUBLIC init_writeField, writefield, close_files CONTAINS SUBROUTINE Writefield_mpi(name_in,field,nind) USE netcdf_mod 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,l,k,n,ncell,q 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 INTEGER :: displ 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_mpi(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 displ=FieldVarId(index)%displ 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)%nc_id(1),Field_val2d, & start=(/ displ+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)%nc_id(1),Field_val3d, & start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /)) DEALLOCATE(field_val3d) ELSE IF (field(1)%ndim==4) THEN DO q=1,FieldVarId(index)%size ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,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)%rval4d(k,:,q) ENDIF ENDDO ENDDO status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d(:,l), & start=(/ displ+ncell,l,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /)) DEALLOCATE(field_val3d) ENDDO ENDIF 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 displ=FieldVarId(index)%displ 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)%nc_id(1), & Field_val2d,start=(/ displ+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)%nc_id(1),Field_val3d, & start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /)) DEALLOCATE(field_val3d) ELSE IF (field(1)%ndim==4) THEN DO q=1,FieldVarId(index)%size ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,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)%rval4d(k+z_down,:,q) 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)%rval4d(k+z_up,:,q) ENDDO ENDDO status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d, & start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /)) DEALLOCATE(field_val3d) ENDDO ENDIF ncell=ncell+n ENDDO ENDIF status=NF90_SYNC(FieldId(Index)) END SUBROUTINE Writefield_mpi SUBROUTINE Create_header_mpi(name,field,nind) USE netcdf_mod USE field_mod USE domain_mod USE spherical_geom_mod USE dimensions USE geometry USE mpi_mod USE mpipara IMPLICIT NONE CHARACTER(LEN=*) :: name CHARACTER(LEN=LEN_TRIM(ADJUSTL(name))) :: name_adj 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,q INTEGER :: iie,jje,iin,jjn INTEGER :: ind_b,ind_e INTEGER :: halo_size LOGICAL :: single INTEGER :: nij INTEGER :: ncell_glo(0:mpi_size-1) INTEGER :: displ, ncell_tot NbField=NbField+1 name_adj=TRIM(ADJUSTL(name)) ! work around ICE with pgf90 FieldName(NbField)=name_adj 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 CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr) displ=0 DO i=1,mpi_rank displ=displ+ncell_glo(i-1) ENDDO FieldVarId(NbField)%displ=displ ncell_tot=sum(ncell_glo(:)) ! status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc', IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid) FieldId(NbField)=ncid status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId) status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) IF (Field(ind_b)%ndim==2) THEN FieldVarId(NbField)%size=1 ALLOCATE(FieldVarId(NbField)%nc_id(1)) ELSE IF (Field(ind_b)%ndim==3) THEN FieldVarId(NbField)%size=1 ALLOCATE(FieldVarId(NbField)%nc_id(1)) status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id) ELSE IF (Field(1)%ndim==4) THEN FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) status = NF90_DEF_DIM(ncid,'lev',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)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/)) ELSE IF (Field(ind_b)%ndim==3) THEN status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, & (/ncell_tot,size(field(ind_b)%rval3d,2),1/)) ELSE IF (Field(ind_b)%ndim==4) THEN DO i=1,FieldVarId(NbField)%size status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /), & FieldVarId(NbField)%nc_id(i)) status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat") status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED, & (/ncell_tot,size(field(ind_b)%rval4d,2),1/)) ENDDO ENDIF 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=(/ displ+ncell /),count=(/ n /)) status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /)) status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+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 CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr) displ=0 DO i=1,mpi_rank displ=displ+ncell_glo(i-1) ENDDO FieldVarId(NbField)%displ=displ ncell_tot=sum(ncell_glo(:)) ! status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc',IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid) ! status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) FieldId(NbField)=ncid status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId) status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) IF (Field(ind_b)%ndim==2) THEN FieldVarId(NbField)%size=1 ALLOCATE(FieldVarId(NbField)%nc_id(1)) ELSE IF (Field(ind_b)%ndim==3) THEN FieldVarId(NbField)%size=1 ALLOCATE(FieldVarId(NbField)%nc_id(1)) status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id) ELSE IF (Field(1)%ndim==4) THEN FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id) 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)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/)) ELSE IF (Field(ind_b)%ndim==3) THEN status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, & (/ncell_tot,size(field(ind_b)%rval3d,2),1/)) ELSE IF (Field(ind_b)%ndim==4) THEN DO q=1,FieldVarId(NbField)%size status = NF90_DEF_VAR(ncid,name_adj//int2str(q),ncprec, & (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q)) status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat") status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED, & (/ncell_tot,size(field(ind_b)%rval4d,2),1/)) ENDDO ENDIF 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 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 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 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 ENDDO ENDDO ENDDO status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /)) status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /)) status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) ncell=ncell+n DEALLOCATE(lon,lat,bounds_lon,bounds_lat) END DO ENDIF END SUBROUTINE Create_Header_mpi end module write_field_mpi_mod