Ignore:
Timestamp:
06/04/19 18:30:08 (5 years ago)
Author:
dubos
Message:

devel : store cell bounds once, use them for XIOS later

File:
1 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/devel/src/output/write_field.f90

    r846 r880  
    11module write_field_mod 
    22  USE genmod 
     3  USE write_field_vars_mod 
    34  IMPLICIT NONE 
    45  PRIVATE   
    5   INTEGER,SAVE :: ncprec 
    6    
    7   TYPE ncvar 
    8     INTEGER :: size 
    9     INTEGER,POINTER :: nc_id(:) 
    10     INTEGER :: displ 
    11   END TYPE ncvar 
    12  
    13   INTEGER, PARAMETER :: MaxWriteField = 1000 
    14   INTEGER, DIMENSION(MaxWriteField),SAVE :: FieldId 
    15   TYPE(ncvar), dimension(MaxWriteField),SAVE :: FieldVarId 
    16   INTEGER, DIMENSION(MaxWriteField),SAVE :: FieldIndex 
    17   CHARACTER(len=255), DIMENSION(MaxWriteField) ::  FieldName  
    18     
    19   INTEGER,SAVE :: NbField = 0 
    206   
    217  PUBLIC init_writeField, writefield, close_files 
    … …  
    4329    END SUBROUTINE init_writeField 
    4430     
    45     function GetFieldIndex(name) 
    46     implicit none 
    47       integer          :: GetFieldindex 
    48       character(len=*) :: name 
    49      
    50       character(len=255) :: TrueName 
    51       integer            :: i 
    52         
    53        
    54       TrueName=TRIM(ADJUSTL(name)) 
    55      
    56       GetFieldIndex=-1 
    57       do i=1,NbField 
    58         if (TrueName==FieldName(i)) then 
    59           GetFieldIndex=i 
    60           exit 
    61         endif 
    62       enddo 
    63     end function GetFieldIndex 
    64      
    65     SUBROUTINE Writefield(name_in,field,nind,once) 
     31   SUBROUTINE Writefield(name_in,field,nind,once) 
    6632    USE domain_mod 
    6733    USE field_mod 
    … …  
    7137    USE netcdf_mod 
    7238    USE grid_param 
    73     IMPLICIT NONE   
    7439     CHARACTER(LEN=*),INTENT(IN) :: name_in 
    7540      TYPE(t_field),POINTER :: field(:) 
    … …  
    8146       
    8247      IF (no_io) RETURN 
    83       IF (grid_type==grid_unst) RETURN 
    8448 
    8549!$OMP BARRIER 
    … …  
    9155      END IF 
    9256 
    93       IF (field(1)%ndim==2) THEN 
    94         CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type) 
    95       ELSE IF (field(1)%ndim==3) THEN 
    96         CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3)  
    97       ELSE IF (field(1)%ndim==4) THEN 
    98         CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3,field(1)%dim4) 
    99       ENDIF 
    100        
    101       CALL gather_field(field,field_glo) 
    102        
     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 
    10371      IF (mpi_rank==0) THEN 
    10472        IF (PRESENT(nind)) THEN 
    … …  
    10977      ENDIF 
    11078       
    111       CALL deallocate_field_glo(field_glo) 
     79      IF(grid_type == grid_ico) CALL deallocate_field_glo(field_glo) 
    11280!$OMP END MASTER 
    11381!$OMP BARRIER 
    … …  
    11583   END SUBROUTINE writefield 
    11684    
    117 !   SUBROUTINE Writefield(name_in,field,nind) 
    118 !   USE netcdf 
    119 !   USE domain_mod 
    120 !   use field_mod 
    121 !   USE dimensions 
    122 !   USE geometry 
    123 !   IMPLICIT NONE 
    124 !     CHARACTER(LEN=*),INTENT(IN) :: name_in 
    125 !     TYPE(t_field),POINTER :: field(:) 
    126 !     INTEGER,OPTIONAL,INTENT(IN) :: nind 
    127 !     REAL(r8),ALLOCATABLE :: field_val2d(:) 
    128 !     REAL(r8),ALLOCATABLE :: field_val3d(:,:) 
    129 !     REAL(r8),ALLOCATABLE :: field_val4d(:,:,:) 
    130 !     TYPE(t_domain),POINTER :: d 
    131 !     INTEGER :: Index 
    132 !     INTEGER :: ind,i,j,k,n,ncell,q 
    133 !     INTEGER :: iie,jje,iin,jjn 
    134 !     INTEGER :: status 
    135 !     CHARACTER(len=255) :: name 
    136 !     CHARACTER(len=255) :: str_ind 
    137 !     INTEGER :: ind_b,ind_e 
    138 !     INTEGER :: halo_size 
    139 !     LOGICAL :: single 
    140 !      
    141 !      
    142 !     name=TRIM(ADJUSTL(name_in)) 
    143  
    144 !     IF (PRESENT(nind)) THEN 
    145 !       name=TRIM(name)//"_"//TRIM(int2str(nind)) 
    146 !       PRINT *,"NAME",nind,int2str(nind),name 
    147 !       ind_b=nind 
    148 !       ind_e=nind 
    149 !       halo_size=1 
    150 !       single=.TRUE. 
    151 !     ELSE 
    152 !       ind_b=1 
    153 !       ind_e=ndomain 
    154 !       halo_size=0 
    155 !       single=.FALSE. 
    156 !     ENDIF       
    157  
    158 !     Index=GetFieldIndex(name) 
    159 !     if (Index==-1) then 
    160 !       call create_header(name,field,nind) 
    161 !       Index=GetFieldIndex(name) 
    162 !     else 
    163 !       FieldIndex(Index)=FieldIndex(Index)+1. 
    164 !     endif 
    165 !      
    166 !     IF (Field(ind_b)%field_type==field_T) THEN 
    167 !       ncell=1 
    168 !       DO ind=ind_b,ind_e 
    169 !         d=>domain(ind) 
    170 !         IF (Field(ind)%field_type/=field_T) THEN 
    171 !           PRINT *,"Writefield, grille non geree" 
    172 !           RETURN 
    173 !         ENDIF 
    174  
    175 !       n=0 
    176 !       DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    177 !         DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    178 !           IF (d%own(i,j) .OR. single) n=n+1 
    179 !         ENDDO 
    180 !       ENDDO 
    181  
    182 !       IF (field(ind)%ndim==2) THEN 
    183 !         ALLOCATE(Field_val2d(n)) 
    184 !         n=0 
    185 !         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    186 !           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    187 !             k=d%iim*(j-1)+i 
    188 !             IF (d%own(i,j) .OR. single) THEN 
    189 !               n=n+1 
    190 !               Field_val2d(n)=field(ind)%rval2d(k) 
    191 !             ENDIF 
    192 !           ENDDO 
    193 !         ENDDO 
    194 !         status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  & 
    195 !                             start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
    196 !         DEALLOCATE(field_val2d) 
    197 !       ELSE IF (field(ind)%ndim==3) THEN 
    198 !         ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2))) 
    199 !         n=0 
    200 !         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    201 !           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    202 !             k=d%iim*(j-1)+i 
    203 !             IF (d%own(i,j) .OR. single) THEN 
    204 !               n=n+1 
    205 !               Field_val3d(n,:)=field(ind)%rval3d(k,:) 
    206 !             ENDIF 
    207 !           ENDDO 
    208 !         ENDDO 
    209 !          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
    210 !                              count=(/n,size(field(1)%rval3d,2),1 /)) 
    211 !         DEALLOCATE(field_val3d) 
    212 !       ELSE IF (field(1)%ndim==4) THEN 
    213  
    214 !         DO q=1,FieldVarId(index)%size 
    215 !          
    216 !           ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2))) 
    217 !           n=0 
    218 !           DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    219 !             DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    220 !               k=d%iim*(j-1)+i 
    221 !               IF (d%own(i,j) .OR. single) THEN 
    222 !                 n=n+1 
    223 !                 Field_val3d(n,:)=field(ind)%rval4d(k,:,q) 
    224 !               ENDIF 
    225 !             ENDDO 
    226 !           ENDDO 
    227 !           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
    228 !                              count=(/n,size(field(1)%rval4d,2),1 /)) 
    229 !           DEALLOCATE(field_val3d) 
    230 !         ENDDO          
    231 !       ENDIF 
    232 !        
    233 !        PRINT *,NF90_STRERROR(status) 
    234 !       ncell=ncell+n 
    235  
    236 !    ENDDO 
    237 !     
    238 !    ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
    239 !       ncell=1 
    240 !       n=0 
    241 !       DO ind=ind_b,ind_e 
    242 !         d=>domain(ind) 
    243 !         CALL swap_geometry(ind) 
    244 !         CALL swap_dimensions(ind) 
    245 !  
    246 !         n=0 
    247 !         DO j=jj_begin+1,jj_end 
    248 !           DO i=ii_begin,ii_end-1 
    249 !             n=n+1 
    250 !           ENDDO 
    251 !         ENDDO 
    252  
    253 !         DO j=jj_begin,jj_end-1 
    254 !           DO i=ii_begin+1,ii_end 
    255 !             n=n+1 
    256 !           ENDDO 
    257 !         ENDDO 
    258  
    259 !       IF (field(ind)%ndim==2) THEN 
    260 !         ALLOCATE(Field_val2d(n)) 
    261  
    262 !         n=0 
    263 !         DO j=jj_begin+1,jj_end 
    264 !           DO i=ii_begin,ii_end-1 
    265 !             n=n+1 
    266 !             k=iim*(j-1)+i 
    267 !             Field_val2d(n)=field(ind)%rval2d(k+z_down) 
    268 !           ENDDO 
    269 !         ENDDO 
    270  
    271 !         DO j=jj_begin,jj_end-1 
    272 !           DO i=ii_begin+1,ii_end 
    273 !             n=n+1 
    274 !             k=iim*(j-1)+i 
    275 !             Field_val2d(n)=field(ind)%rval2d(k+z_up) 
    276 !           ENDDO 
    277 !         ENDDO 
    278  
    279 !         status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       & 
    280 !                             Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
    281 !         DEALLOCATE(field_val2d) 
    282  
    283 !       ELSE IF (field(ind)%ndim==3) THEN 
    284 !         ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2))) 
    285 !         n=0 
    286 !         DO j=jj_begin+1,jj_end 
    287 !           DO i=ii_begin,ii_end-1 
    288 !             n=n+1 
    289 !             k=iim*(j-1)+i 
    290 !             Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:) 
    291 !           ENDDO 
    292 !         ENDDO 
    293  
    294 !         DO j=jj_begin,jj_end-1 
    295 !           DO i=ii_begin+1,ii_end 
    296 !             n=n+1 
    297 !             k=iim*(j-1)+i 
    298 !             Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:) 
    299 !           ENDDO 
    300 !         ENDDO 
    301 !          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
    302 !                              count=(/n,size(field(1)%rval3d,2),1 /)) 
    303 !         DEALLOCATE(field_val3d) 
    304 !       ELSE IF (field(1)%ndim==4) THEN 
    305  
    306 !         DO q=1,FieldVarId(index)%size 
    307 !           ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2))) 
    308 !           n=0 
    309 !           DO j=jj_begin+1,jj_end 
    310 !             DO i=ii_begin,ii_end-1 
    311 !               n=n+1 
    312 !               k=iim*(j-1)+i 
    313 !               Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q) 
    314 !             ENDDO 
    315 !           ENDDO 
    316  
    317 !           DO j=jj_begin,jj_end-1 
    318 !             DO i=ii_begin+1,ii_end 
    319 !               n=n+1 
    320 !               k=iim*(j-1)+i 
    321 !               Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q) 
    322 !             ENDDO 
    323 !           ENDDO 
    324  
    325 !           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,1,FieldIndex(Index) /),   & 
    326 !                               count=(/n,size(field(1)%rval4d,2),1 /)) 
    327 !           DEALLOCATE(field_val3d) 
    328 !         ENDDO 
    329 !       ENDIF 
    330 !        
    331 !        PRINT *,NF90_STRERROR(status) 
    332 !       ncell=ncell+n 
    333  
    334 !    ENDDO 
    335 !     
    336 !    ENDIF 
    337 !    status=NF90_SYNC(FieldId(Index)) 
    338 !      
    339 !   END SUBROUTINE Writefield 
    340  
    341  
    342     SUBROUTINE Writefield_gen(name_in, field, domain_type, ind_b_in, ind_e_in,once ) 
    343     USE netcdf_mod 
    344     USE domain_mod 
    345     USE field_mod 
    346     IMPLICIT NONE 
    347       CHARACTER(LEN=*),INTENT(IN) :: name_in 
    348       TYPE(t_field), POINTER :: field(:) 
    349       TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:) 
    350       INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in 
    351       INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in 
    352       REAL(r8),ALLOCATABLE :: field_val2d(:) 
    353       REAL(r8),ALLOCATABLE :: field_val3d(:,:) 
    354       REAL(r8),ALLOCATABLE :: field_val4d(:,:,:) 
    355       LOGICAL, INTENT(IN) :: once 
    356       TYPE(t_domain),POINTER :: d 
    357       INTEGER :: Index 
    358       INTEGER :: ind,i,j,k,n,ncell,q 
    359       INTEGER :: iie,jje,iin,jjn 
    360       INTEGER :: status 
    361       CHARACTER(len=255) :: name 
    362       CHARACTER(len=255) :: str_ind 
    363       INTEGER :: ind_b,ind_e 
    364       INTEGER :: halo_size 
    365       LOGICAL :: single 
    366        
    367       name=TRIM(ADJUSTL(name_in)) 
     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)) 
    368108 
    369109      IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN 
    … …  
    371111        ind_b=ind_b_in 
    372112        ind_e=ind_b_in 
    373         halo_size=1 
    374113        single=.TRUE. 
    375114      ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
    … …  
    377116        ind_b=ind_e_in 
    378117        ind_e=ind_e_in 
    379         halo_size=1 
    380118        single=.TRUE. 
    381119      ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
    382120        ind_b=ind_b_in 
    383121        ind_e=ind_e_in 
    384         halo_size=0 
    385122        single=.FALSE. 
    386123      ELSE 
    387         halo_size=0 
    388124        ind_b=1 
    389125        ind_e=ndomain 
    390126        single=.FALSE. 
    391127      ENDIF       
    392        
    393       Index=GetFieldIndex(name) 
    394       if (Index==-1) then 
    395         call create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once) 
    396         Index=GetFieldIndex(name) 
     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) 
    397133      else 
    398         FieldIndex(Index)=FieldIndex(Index)+1. 
     134         FieldIndex(index)=FieldIndex(index)+1. 
    399135      endif 
    400        
    401       IF (Field(ind_b)%field_type==field_T) THEN 
    402  
    403         ncell=0 
    404         DO ind=ind_b,ind_e 
    405           d=>domain_type(ind) 
    406           DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    407             DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    408               IF (d%assign_domain(i,j)==ind .OR. single) ncell=ncell+1 
    409             ENDDO 
    410           ENDDO 
    411         ENDDO 
    412          
    413         IF (field(1)%ndim==2) THEN 
    414           ALLOCATE(Field_val2d(ncell)) 
    415           n=0 
     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             PRINT *, 'write_2d :', ind, n, n_begin+n, ij  
     178             field_val2d(n_begin+n) = field(ind)%rval2d(ij) 
     179          END DO 
     180          n_begin = n_begin + cells(ind)%ncell 
     181       END DO 
     182       IF (once) THEN  
     183          status=NF90_PUT_VAR(FieldId(index), FieldVarId(index)%nc_id(1), & 
     184               Field_val2d, start=(/ 1 /),count=(/ncell /)) 
     185       ELSE 
     186          status=NF90_PUT_VAR(FieldId(Index), FieldVarId(index)%nc_id(1), & 
     187               Field_val2d,  start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /)) 
     188       ENDIF 
     189       DEALLOCATE(field_val2d) 
     190     END SUBROUTINE write_2d 
     191 
     192     SUBROUTINE write_3d 
     193       dim3 = SIZE(field(1)%rval3d,2) 
     194       ALLOCATE(field_val3d(ncell,dim3)) 
     195       n_begin=0   
     196       DO ind=ind_b,ind_e 
     197          DO n=1, cells(ind)%ncell 
     198             ij=cells(ind)%ij(n) 
     199             field_val3d(n_begin+n,:) = field(ind)%rval3d(ij,:) 
     200          END DO 
     201       END DO 
     202       IF (once) THEN  
     203          status=NF90_PUT_VAR(FieldId(Index), FieldVarId(index)%nc_id(1), & 
     204               field_val3d, start=(/ 1,1 /), & 
     205               count=(/ncell,dim3 /)) 
     206       ELSE 
     207          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1), & 
     208               field_val3d, start=(/ 1,1,FieldIndex(Index) /),   & 
     209               count=(/ncell,size(field(1)%rval3d,2),1 /)) 
     210       ENDIF 
     211       DEALLOCATE(field_val3d) 
     212     END SUBROUTINE write_3d 
     213 
     214     SUBROUTINE write_4d 
     215       dim3 = SIZE(field(1)%rval4d,2) 
     216       ALLOCATE(field_val3d(ncell,dim3)) 
     217       DO q=1,FieldVarId(index)%size 
     218          n_begin=0   
    416219          DO ind=ind_b,ind_e 
    417             d=>domain_type(ind) 
    418             DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    419               DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    420                 k=d%iim*(j-1)+i 
    421                 IF (d%assign_domain(i,j)==ind .OR. single) THEN 
    422                   n=n+1 
    423                   Field_val2d(n)=field(ind)%rval2d(k) 
    424                 ENDIF 
    425               ENDDO 
    426             ENDDO 
    427           ENDDO 
    428  
     220             DO n=1, cells(ind)%ncell 
     221                ij=cells(ind)%ij(n) 
     222                field_val3d(n_begin+n,:) = field(ind)%rval4d(ij,:,q) 
     223             END DO 
     224          END DO 
    429225          IF (once) THEN  
    430             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  & 
    431                                start=(/ 1 /),count=(/ncell /)) 
     226             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q), & 
     227                  field_val3d, start=(/ 1,1 /),   & 
     228                  count=(/ncell,dim3,1 /)) 
    432229          ELSE 
    433             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  & 
    434                                start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /)) 
     230             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q), & 
     231                  field_val3d, start=(/ 1,1,FieldIndex(Index) /),   & 
     232                  count=(/ncell,dim3,1 /)) 
    435233          ENDIF 
    436            
    437           DEALLOCATE(field_val2d) 
    438         ELSE IF (field(1)%ndim==3) THEN 
    439           ALLOCATE(Field_val3d(ncell,size(field(1)%rval3d,2))) 
    440           n=0 
    441           DO ind=ind_b,ind_e 
    442             d=>domain_type(ind) 
    443             DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    444               DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    445                 k=d%iim*(j-1)+i 
    446                 IF (d%assign_domain(i,j)==ind .OR. single) THEN 
    447                   n=n+1 
    448                   Field_val3d(n,:)=field(ind)%rval3d(k,:) 
    449                 ENDIF 
    450               ENDDO 
    451             ENDDO 
    452           ENDDO 
    453  
    454           PRINT *, 'Writefield ', TRIM(name), MAXVAL(Field_val3d(:,1)), MINVAL(Field_val3d(:,1)) ! FIXME 
    455  
    456           IF (once) THEN  
    457             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1 /),   & 
    458                                  count=(/ncell,size(field(1)%rval3d,2) /)) 
    459           ELSE 
    460              status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1,FieldIndex(Index) /),   & 
    461                                  count=(/ncell,size(field(1)%rval3d,2),1 /)) 
    462           ENDIF 
    463                                  
    464           DEALLOCATE(field_val3d) 
    465         ELSE IF (field(1)%ndim==4) THEN 
    466  
    467           DO q=1,FieldVarId(index)%size 
    468            
    469             ALLOCATE(Field_val3d(ncell,size(field(1)%rval4d,2))) 
    470             n=0 
    471             DO ind=ind_b,ind_e 
    472               d=>domain_type(ind) 
    473               DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    474                 DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    475                   k=d%iim*(j-1)+i 
    476                   IF (d%assign_domain(i,j)==ind .OR. single) THEN 
    477                     n=n+1 
    478                     Field_val3d(n,:)=field(ind)%rval4d(k,:,q) 
    479                   ENDIF 
    480                 ENDDO 
    481               ENDDO 
    482             ENDDO 
    483              
    484             IF (once) THEN  
    485               status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1 /),   & 
    486                                  count=(/ncell,size(field(1)%rval4d,2),1 /)) 
    487             ELSE 
    488               status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,FieldIndex(Index) /),   & 
    489                                  count=(/ncell,size(field(1)%rval4d,2),1 /)) 
    490             ENDIF 
    491             DEALLOCATE(field_val3d) 
    492           ENDDO          
    493         ENDIF 
    494          
    495       
    496      ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
    497  
    498         ncell=0 
    499         DO ind=ind_b,ind_e 
    500           d=>domain_type(ind) 
    501    
    502           DO j=d%jj_begin+1,d%jj_end 
    503             DO i=d%ii_begin,d%ii_end-1 
    504               ncell=ncell+1 
    505             ENDDO 
    506           ENDDO 
    507  
    508           DO j=d%jj_begin,d%jj_end-1 
    509             DO i=d%ii_begin+1,d%ii_end 
    510               ncell=ncell+1 
    511             ENDDO 
    512           ENDDO 
    513         ENDDO 
    514  
    515         IF (field(1)%ndim==2) THEN 
    516           ALLOCATE(Field_val2d(ncell)) 
    517  
    518           n=0 
    519            
    520           DO ind=ind_b,ind_e 
    521             d=>domain_type(ind) 
    522             DO j=d%jj_begin+1,d%jj_end 
    523               DO i=d%ii_begin,d%ii_end-1 
    524                 n=n+1 
    525                 k=d%iim*(j-1)+i 
    526                 Field_val2d(n)=field(ind)%rval2d(k+d%z_down) 
    527               ENDDO 
    528             ENDDO 
    529  
    530             DO j=d%jj_begin,d%jj_end-1 
    531               DO i=d%ii_begin+1,d%ii_end 
    532                 n=n+1 
    533                 k=d%iim*(j-1)+i 
    534                 Field_val2d(n)=field(ind)%rval2d(k+d%z_up) 
    535               ENDDO 
    536             ENDDO 
    537           ENDDO 
    538            
    539           IF (once) THEN  
    540             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       & 
    541                                 Field_val2d,start=(/ 1 /),count=(/ncell /)) 
    542            ELSE 
    543             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       & 
    544                                 Field_val2d,start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /)) 
    545           ENDIF 
    546           DEALLOCATE(field_val2d) 
    547  
    548         ELSE IF (field(1)%ndim==3) THEN 
    549           ALLOCATE(Field_val3d(ncell,size(field(1)%rval3d,2))) 
    550           n=0 
    551           DO ind=ind_b,ind_e 
    552             d=>domain_type(ind) 
    553             DO j=d%jj_begin+1,d%jj_end 
    554               DO i=d%ii_begin,d%ii_end-1 
    555                 n=n+1 
    556                 k=d%iim*(j-1)+i 
    557                 Field_val3d(n,:)=field(ind)%rval3d(k+d%z_down,:) 
    558               ENDDO 
    559             ENDDO 
    560  
    561             DO j=d%jj_begin,d%jj_end-1 
    562               DO i=d%ii_begin+1,d%ii_end 
    563                 n=n+1 
    564                 k=d%iim*(j-1)+i 
    565                 Field_val3d(n,:)=field(ind)%rval3d(k+d%z_up,:) 
    566               ENDDO 
    567             ENDDO 
    568           ENDDO 
    569            
    570           IF (once) THEN  
    571             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1 /),   & 
    572                                 count=(/ncell,size(field(1)%rval3d,2) /)) 
    573           ELSE 
    574             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1,FieldIndex(Index) /),   & 
    575                                 count=(/ncell,size(field(1)%rval3d,2),1 /)) 
    576           ENDIF 
    577           DEALLOCATE(field_val3d) 
    578          
    579         ELSE IF (field(1)%ndim==4) THEN 
    580  
    581           DO q=1,FieldVarId(index)%size 
    582             ALLOCATE(Field_val3d(ncell,size(field(1)%rval4d,2))) 
    583             n=0 
    584             DO ind=ind_b,ind_e 
    585               d=>domain_type(ind) 
    586               DO j=d%jj_begin+1,d%jj_end 
    587                 DO i=d%ii_begin,d%ii_end-1 
    588                   n=n+1 
    589                   k=d%iim*(j-1)+i 
    590                   Field_val3d(n,:)=field(ind)%rval4d(k+d%z_down,:,q) 
    591                 ENDDO 
    592               ENDDO 
    593  
    594               DO j=d%jj_begin,d%jj_end-1 
    595                 DO i=d%ii_begin+1,d%ii_end 
    596                   n=n+1 
    597                   k=d%iim*(j-1)+i 
    598                   Field_val3d(n,:)=field(ind)%rval4d(k+d%z_up,:,q) 
    599                 ENDDO 
    600               ENDDO 
    601             ENDDO 
    602  
    603             IF (once) THEN  
    604               status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,1 /),   & 
    605                                   count=(/ncell,size(field(1)%rval4d,2) /)) 
    606             ELSE 
    607               status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,1,FieldIndex(Index) /),   & 
    608                                   count=(/ncell,size(field(1)%rval4d,2),1 /)) 
    609             ENDIF 
    610             DEALLOCATE(field_val3d) 
    611           ENDDO 
    612         ENDIF 
    613       
    614      ENDIF 
    615      status=NF90_SYNC(FieldId(Index)) 
    616        
    617     END SUBROUTINE Writefield_gen 
    618  
    619        
    620  
    621     SUBROUTINE Writefield_mpi(name_in,field,nind) 
    622     USE netcdf_mod 
    623     USE domain_mod 
    624     use field_mod 
    625     USE dimensions 
    626     USE geometry 
    627     IMPLICIT NONE 
    628       CHARACTER(LEN=*),INTENT(IN) :: name_in 
    629       TYPE(t_field),POINTER :: field(:) 
    630       INTEGER,OPTIONAL,INTENT(IN) :: nind 
    631       REAL(r8),ALLOCATABLE :: field_val2d(:) 
    632       REAL(r8),ALLOCATABLE :: field_val3d(:,:) 
    633       REAL(r8),ALLOCATABLE :: field_val4d(:,:,:) 
    634       TYPE(t_domain),POINTER :: d 
    635       INTEGER :: Index 
    636       INTEGER :: ind,i,j,l,k,n,ncell,q 
    637       INTEGER :: iie,jje,iin,jjn 
    638       INTEGER :: status 
    639       CHARACTER(len=255) :: name 
    640       CHARACTER(len=255) :: str_ind 
    641       INTEGER :: ind_b,ind_e 
    642       INTEGER :: halo_size 
    643       LOGICAL :: single 
    644       INTEGER :: displ 
    645        
    646        
    647       name=TRIM(ADJUSTL(name_in)) 
    648  
    649       IF (PRESENT(nind)) THEN 
    650         name=TRIM(name)//"_"//TRIM(int2str(nind)) 
    651         PRINT *,"NAME",nind,int2str(nind),name 
    652         ind_b=nind 
    653         ind_e=nind 
    654         halo_size=1 
    655         single=.TRUE. 
    656       ELSE 
    657         ind_b=1 
    658         ind_e=ndomain 
    659         halo_size=0 
    660         single=.FALSE. 
    661       ENDIF       
    662  
    663       Index=GetFieldIndex(name) 
    664       if (Index==-1) then 
    665         call create_header_mpi(name,field,nind) 
    666         Index=GetFieldIndex(name) 
    667       else 
    668         FieldIndex(Index)=FieldIndex(Index)+1. 
    669       endif 
    670        
    671       IF (Field(ind_b)%field_type==field_T) THEN 
    672         ncell=1 
    673         DO ind=ind_b,ind_e 
    674           d=>domain(ind) 
    675           IF (Field(ind)%field_type/=field_T) THEN 
    676             PRINT *,"Writefield, grille non geree" 
    677             RETURN 
    678           ENDIF 
    679  
    680         n=0 
    681         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    682           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    683             IF (d%own(i,j) .OR. single) n=n+1 
    684           ENDDO 
    685         ENDDO 
    686  
    687         displ=FieldVarId(index)%displ 
    688          
    689         IF (field(ind)%ndim==2) THEN 
    690           ALLOCATE(Field_val2d(n)) 
    691           n=0 
    692           DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    693             DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    694               k=d%iim*(j-1)+i 
    695               IF (d%own(i,j) .OR. single) THEN 
    696                 n=n+1 
    697                 Field_val2d(n)=field(ind)%rval2d(k) 
    698               ENDIF 
    699             ENDDO 
    700           ENDDO 
    701           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  & 
    702                               start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /)) 
    703           DEALLOCATE(field_val2d) 
    704         ELSE IF (field(ind)%ndim==3) THEN 
    705           ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2))) 
    706           n=0 
    707           DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    708             DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    709               k=d%iim*(j-1)+i 
    710               IF (d%own(i,j) .OR. single) THEN 
    711                 n=n+1 
    712                 Field_val3d(n,:)=field(ind)%rval3d(k,:) 
    713               ENDIF 
    714             ENDDO 
    715           ENDDO 
    716           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,  & 
    717                                 start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /)) 
    718           DEALLOCATE(field_val3d) 
    719  
    720         ELSE IF (field(1)%ndim==4) THEN 
    721  
    722           DO q=1,FieldVarId(index)%size 
    723            
    724             ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2))) 
    725             n=0 
    726             DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    727               DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    728                 k=d%iim*(j-1)+i 
    729                 IF (d%own(i,j) .OR. single) THEN 
    730                   n=n+1 
    731                   Field_val3d(n,:)=field(ind)%rval4d(k,:,q) 
    732                 ENDIF 
    733               ENDDO 
    734              ENDDO 
    735  
    736             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d(:,l), & 
    737                                 start=(/ displ+ncell,l,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /)) 
    738             DEALLOCATE(field_val3d) 
    739           ENDDO          
    740         ENDIF 
    741          
    742         ncell=ncell+n 
    743  
    744      ENDDO 
    745       
    746      ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
    747         ncell=1 
    748         n=0 
    749         DO ind=ind_b,ind_e 
    750           d=>domain(ind) 
    751           CALL swap_geometry(ind) 
    752           CALL swap_dimensions(ind) 
    753    
    754           n=0 
    755           DO j=jj_begin+1,jj_end 
    756             DO i=ii_begin,ii_end-1 
    757               n=n+1 
    758             ENDDO 
    759           ENDDO 
    760  
    761           DO j=jj_begin,jj_end-1 
    762             DO i=ii_begin+1,ii_end 
    763               n=n+1 
    764             ENDDO 
    765           ENDDO 
    766  
    767         displ=FieldVarId(index)%displ 
    768  
    769         IF (field(ind)%ndim==2) THEN 
    770           ALLOCATE(Field_val2d(n)) 
    771  
    772           n=0 
    773           DO j=jj_begin+1,jj_end 
    774             DO i=ii_begin,ii_end-1 
    775               n=n+1 
    776               k=iim*(j-1)+i 
    777               Field_val2d(n)=field(ind)%rval2d(k+z_down) 
    778             ENDDO 
    779           ENDDO 
    780  
    781           DO j=jj_begin,jj_end-1 
    782             DO i=ii_begin+1,ii_end 
    783               n=n+1 
    784               k=iim*(j-1)+i 
    785               Field_val2d(n)=field(ind)%rval2d(k+z_up) 
    786             ENDDO 
    787           ENDDO 
    788  
    789           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       & 
    790                               Field_val2d,start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /)) 
    791           DEALLOCATE(field_val2d) 
    792  
    793         ELSE IF (field(ind)%ndim==3) THEN 
    794           ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2))) 
    795           n=0 
    796           DO j=jj_begin+1,jj_end 
    797             DO i=ii_begin,ii_end-1 
    798               n=n+1 
    799               k=iim*(j-1)+i 
    800               Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:) 
    801             ENDDO 
    802           ENDDO 
    803  
    804           DO j=jj_begin,jj_end-1 
    805             DO i=ii_begin+1,ii_end 
    806               n=n+1 
    807               k=iim*(j-1)+i 
    808               Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:) 
    809             ENDDO 
    810           ENDDO 
    811  
    812            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,    & 
    813                                start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /)) 
    814           DEALLOCATE(field_val3d) 
    815  
    816         ELSE IF (field(1)%ndim==4) THEN 
    817  
    818           DO q=1,FieldVarId(index)%size 
    819             ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2))) 
    820             n=0 
    821             DO j=jj_begin+1,jj_end 
    822               DO i=ii_begin,ii_end-1 
    823                 n=n+1 
    824                 k=iim*(j-1)+i 
    825                 Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q) 
    826               ENDDO 
    827             ENDDO 
    828  
    829             DO j=jj_begin,jj_end-1 
    830               DO i=ii_begin+1,ii_end 
    831                 n=n+1 
    832                 k=iim*(j-1)+i 
    833                 Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q) 
    834               ENDDO 
    835             ENDDO 
    836  
    837             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,  & 
    838                                 start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /)) 
    839             DEALLOCATE(field_val3d) 
    840           ENDDO 
    841         ENDIF 
    842          
    843         ncell=ncell+n 
    844  
    845      ENDDO 
    846       
    847      ENDIF 
    848      status=NF90_SYNC(FieldId(Index)) 
    849        
    850     END SUBROUTINE Writefield_mpi 
    851      
    852      
    853             
    854 !   SUBROUTINE Create_header(name,field,nind) 
    855 !   USE netcdf 
    856 !   USE field_mod 
    857 !   USE domain_mod 
    858 !   USE spherical_geom_mod 
    859 !   USE dimensions 
    860 !   USE geometry 
    861 !   IMPLICIT NONE 
    862 !     CHARACTER(LEN=*) :: name 
    863 !     TYPE(t_field),POINTER :: field(:) 
    864 !     INTEGER,OPTIONAL,INTENT(IN) :: nind 
    865 !     INTEGER :: ncell 
    866 !     INTEGER :: nvert 
    867 !     REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:) 
    868 !     TYPE(t_domain),POINTER :: d 
    869 !     INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId 
    870 !     INTEGER :: dim3id,dim4id 
    871 !     INTEGER :: status 
    872 !     INTEGER :: ind,i,j,k,n,q 
    873 !     INTEGER :: iie,jje,iin,jjn 
    874 !     INTEGER :: ind_b,ind_e 
    875 !     INTEGER :: halo_size 
    876 !     LOGICAL :: single  
    877 !     INTEGER :: nij 
    878 !          
    879 !     NbField=NbField+1 
    880 !     FieldName(NbField)=TRIM(ADJUSTL(name)) 
    881 !     FieldIndex(NbField)=1 
    882 !      
    883 !     IF (PRESENT(nind)) THEN 
    884 !       ind_b=nind 
    885 !       ind_e=nind 
    886 !       halo_size=1 
    887 !       single=.TRUE. 
    888 !     ELSE 
    889 !       ind_b=1 
    890 !       ind_e=ndomain 
    891 !       halo_size=0 
    892 !       single=.FALSE. 
    893 !     ENDIF 
    894 !      
    895 !     ncell=0 
    896 !      
    897 !     IF (Field(ind_b)%field_type==field_T) THEN 
    898 !       nvert=6 
    899 !        
    900 !       DO ind=ind_b,ind_e 
    901 !         d=>domain(ind) 
    902 !        
    903 !         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    904 !           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    905 !             IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1 
    906 !           ENDDO 
    907 !         ENDDO 
    908  
    909 !       END DO 
    910 !      
    911 !       status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
    912 !       FieldId(NbField)=ncid 
    913 !       status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
    914 !       status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
    915  
    916 !       IF (Field(ind_b)%ndim==2)  THEN 
    917 !         FieldVarId(NbField)%size=1 
    918 !         ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    919 !       ELSE IF (Field(ind_b)%ndim==3)  THEN 
    920 !         FieldVarId(NbField)%size=1 
    921 !         ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    922 !         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) 
    923 !       ELSE IF (Field(1)%ndim==4) THEN 
    924 !         FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
    925 !         ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
    926 !         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) 
    927 !          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id) 
    928 !       ENDIF 
    929 !      
    930 !       status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
    931 !      
    932 !       status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) 
    933 !       status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
    934 !       status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
    935 !       status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") 
    936 !       status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) 
    937 !       status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
    938 !       status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
    939 !       status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") 
    940 !       status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
    941 !       status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
    942  
    943 !       IF (Field(ind_b)%ndim==2) THEN 
    944 !         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
    945 !         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
    946 !       ELSE IF (Field(ind_b)%ndim==3) THEN 
    947 !         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
    948 !         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
    949 !       ELSE IF (Field(ind_b)%ndim==4) THEN 
    950 !         DO i=1,FieldVarId(NbField)%size 
    951 !           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),  & 
    952 !                                 FieldVarId(NbField)%nc_id(i)) 
    953 !           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat") 
    954 !         ENDDO         
    955 !       ENDIF  
    956 !  
    957 !      
    958 !       status = NF90_ENDDEF(ncid)       
    959  
    960 !       ncell=1 
    961 !       DO ind=ind_b,ind_e 
    962 !         d=>domain(ind) 
    963 !  
    964 !         n=0 
    965 !         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    966 !           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    967 !             IF (single .OR. d%own(i,j)) n=n+1 
    968 !           ENDDO 
    969 !         ENDDO 
    970 !          
    971 !        ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n)) 
    972 !          
    973 !         n=0   
    974 !         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    975 !           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    976 !               IF (d%own(i,j) .OR. single) THEN 
    977 !               n=n+1 
    978 !               CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n)) 
    979 !               lon(n)=lon(n)*180/Pi 
    980 !               lat(n)=lat(n)*180/Pi 
    981 !               DO k=0,5 
    982 !                 CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n)) 
    983 !                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
    984 !                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
    985 !               ENDDO 
    986 !             ENDIF 
    987 !           ENDDO 
    988 !         ENDDO 
    989 !         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /)) 
    990 !         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /)) 
    991 !         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
    992 !         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
    993 !  
    994 !         ncell=ncell+n 
    995 !         DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
    996 !     END DO  
    997  
    998 !   ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
    999 !       nvert=3 
    1000 !       DO ind=ind_b,ind_e 
    1001 !         d=>domain(ind) 
    1002 !        
    1003 !         DO j=d%jj_begin+1,d%jj_end 
    1004 !           DO i=d%ii_begin,d%ii_end-1 
    1005 !             ncell=ncell+1 
    1006 !           ENDDO 
    1007 !         ENDDO 
    1008  
    1009 !         DO j=d%jj_begin,d%jj_end-1 
    1010 !           DO i=d%ii_begin+1,d%ii_end 
    1011 !             ncell=ncell+1 
    1012 !           ENDDO 
    1013 !         ENDDO 
    1014  
    1015 !       END DO 
    1016 !      
    1017 !       status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
    1018 !       FieldId(NbField)=ncid 
    1019 !       status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
    1020 !       status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
    1021  
    1022 !       IF (Field(ind_b)%ndim==2)  THEN 
    1023 !         FieldVarId(NbField)%size=1 
    1024 !         ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    1025 !       ELSE IF (Field(ind_b)%ndim==3)  THEN 
    1026 !         FieldVarId(NbField)%size=1 
    1027 !         ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    1028 !         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) 
    1029 !       ELSE IF (Field(1)%ndim==4) THEN 
    1030 !         FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
    1031 !         ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
    1032 !         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) 
    1033 !       ENDIF 
    1034  
    1035  
    1036 !      
    1037 !       status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
    1038 !      
    1039 !       status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) 
    1040 !       status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
    1041 !       status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
    1042 !       status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") 
    1043 !       status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) 
    1044 !       status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
    1045 !       status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
    1046 !       status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") 
    1047 !       status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
    1048 !       status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
    1049  
    1050  
    1051 !       IF (Field(ind_b)%ndim==2) THEN 
    1052 !         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
    1053 !         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
    1054 !       ELSE IF (Field(ind_b)%ndim==3) THEN 
    1055 !         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
    1056 !         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
    1057 !       ELSE IF (Field(ind_b)%ndim==4) THEN 
    1058 !         DO q=1,FieldVarId(NbField)%size 
    1059 !           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),NF90_DOUBLE,             & 
    1060 !                                 (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q)) 
    1061 !           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat") 
    1062 !         ENDDO         
    1063 !       ENDIF  
    1064 !        
    1065 !       status = NF90_ENDDEF(ncid)       
    1066  
    1067 !       ncell=1 
    1068 !       DO ind=ind_b,ind_e 
    1069 !         d=>domain(ind) 
    1070 !         CALL swap_geometry(ind) 
    1071 !         CALL swap_dimensions(ind) 
    1072 !  
    1073 !         n=0 
    1074 !         DO j=jj_begin+1,jj_end 
    1075 !           DO i=ii_begin,ii_end-1 
    1076 !             n=n+1 
    1077 !           ENDDO 
    1078 !         ENDDO 
    1079  
    1080 !         DO j=jj_begin,jj_end-1 
    1081 !           DO i=ii_begin+1,ii_end 
    1082 !             n=n+1 
    1083 !           ENDDO 
    1084 !         ENDDO 
    1085  
    1086 !        ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n)) 
    1087 !          
    1088 !         n=0   
    1089 !       
    1090 !         DO j=jj_begin+1,jj_end 
    1091 !           DO i=ii_begin,ii_end-1 
    1092 !             nij=(j-1)*iim+i 
    1093 !             n=n+1 
    1094 !             CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n)) 
    1095 !             lon(n)=lon(n)*180/Pi 
    1096 !!             IF (lon(n)<0) lon(n)=lon(n)+360 
    1097 !             lat(n)=lat(n)*180/Pi 
    1098 !             CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
    1099 !             CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) 
    1100 !             CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) 
    1101  
    1102 !             DO k=0,2 
    1103 !               bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
    1104 !               bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
    1105 !                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 
    1106 !             ENDDO 
    1107 !           ENDDO 
    1108 !         ENDDO 
    1109  
    1110 !         DO j=jj_begin,jj_end-1 
    1111 !           DO i=ii_begin+1,ii_end 
    1112 !             nij=(j-1)*iim+i 
    1113 !             n=n+1 
    1114 !             CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n)) 
    1115 !             lon(n)=lon(n)*180/Pi 
    1116 !              IF (lon(n)<0) lon(n)=lon(n)+360 
    1117 !             lat(n)=lat(n)*180/Pi 
    1118 !             CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
    1119 !             CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) 
    1120 !             CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) 
    1121  
    1122 !             DO k=0,2 
    1123 !               bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
    1124 !               bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
    1125 !                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 
    1126 !             ENDDO 
    1127 !           ENDDO 
    1128 !         ENDDO 
    1129 !          
    1130 !          
    1131 !         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /)) 
    1132 !         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /)) 
    1133 !         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
    1134 !         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
    1135 !  
    1136 !         ncell=ncell+n 
    1137 !         DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
    1138 !     END DO           
    1139 !   ENDIF 
    1140  
    1141  
    1142 !      
    1143 !      status = NF90_CLOSE(ncid) 
    1144  
    1145 !   END SUBROUTINE Create_Header 
    1146  
    1147  
    1148  
     234       END DO 
     235       DEALLOCATE(field_val3d) 
     236     END SUBROUTINE write_4d 
     237 
     238   END SUBROUTINE Writefield_gen 
    1149239    
    1150     SUBROUTINE Create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once) 
    1151     USE netcdf_mod 
    1152     USE field_mod 
    1153     USE domain_mod 
    1154     USE metric 
    1155     USE spherical_geom_mod 
    1156     IMPLICIT NONE 
    1157       CHARACTER(LEN=*),INTENT(IN) :: name_in 
    1158       TYPE(t_field),POINTER :: field(:) 
    1159       TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:) 
    1160       INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in 
    1161       INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in 
    1162       LOGICAL,INTENT(IN) :: once 
    1163        
    1164       INTEGER :: ncell 
    1165       INTEGER :: nvert 
    1166       REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:) 
    1167       TYPE(t_domain),POINTER :: d 
    1168       INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId 
    1169       INTEGER :: dim3id,dim4id 
    1170       INTEGER :: status 
    1171       INTEGER :: ind,i,j,k,n,q 
    1172       INTEGER :: iie,jje,iin,jjn 
    1173       INTEGER :: ind_b,ind_e 
    1174       INTEGER :: halo_size 
    1175       LOGICAL :: single  
    1176       INTEGER :: nij 
    1177       CHARACTER(LEN=255) :: name 
    1178       INTEGER :: l,level_size, levId, dimlevId 
    1179              
    1180       name=TRIM(ADJUSTL(name_in)) 
    1181  
    1182       IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN 
     240   SUBROUTINE Create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once) 
     241     USE netcdf_mod 
     242     USE field_mod 
     243     USE domain_mod 
     244     USE metric 
     245     USE spherical_geom_mod 
     246     CHARACTER(LEN=*),INTENT(IN) :: name_in 
     247     TYPE(t_field),POINTER :: field(:) 
     248     TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:) 
     249     INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in 
     250     INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in 
     251     LOGICAL,INTENT(IN) :: once 
     252      
     253     TYPE(t_cellset), POINTER :: cells(:) 
     254     INTEGER :: ncell 
     255     INTEGER :: nvert 
     256     REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:) 
     257     TYPE(t_domain),POINTER :: d 
     258     INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId 
     259     INTEGER :: dim3id,dim4id 
     260     INTEGER :: status 
     261     INTEGER :: ind,i,j,k,n,q 
     262     INTEGER :: iie,jje,iin,jjn 
     263     INTEGER :: ind_b,ind_e, n_begin, n_end 
     264     LOGICAL :: single  
     265     INTEGER :: nij 
     266     CHARACTER(LEN=255) :: name 
     267     INTEGER :: l,level_size, levId, dimlevId 
     268      
     269     name=TRIM(ADJUSTL(name_in)) 
     270      
     271     IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN 
    1183272        name=TRIM(name)//"_"//TRIM(int2str(ind_b)) 
    1184273        ind_b=ind_b_in 
    1185274        ind_e=ind_b_in 
    1186         halo_size=1 
    1187275        single=.TRUE. 
    1188       ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
     276     ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
    1189277        name=TRIM(name)//"_"//TRIM(int2str(ind_e)) 
    1190278        ind_b=ind_e_in 
    1191279        ind_e=ind_e_in 
    1192         halo_size=1 
    1193280        single=.TRUE. 
    1194       ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
     281     ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
    1195282        ind_b=ind_b_in 
    1196283        ind_e=ind_e_in 
    1197         halo_size=0 
    1198284        single=.FALSE. 
    1199       ELSE 
    1200         halo_size=0 
     285     ELSE 
    1201286        ind_b=1 
    1202287        ind_e=ndomain 
    1203288        single=.FALSE. 
    1204       ENDIF       
    1205                  
    1206       NbField=NbField+1 
    1207       FieldName(NbField)=TRIM(ADJUSTL(name)) 
    1208       FieldIndex(NbField)=1 
    1209        
    1210       ncell=0 
    1211        
    1212       IF (Field(ind_b)%field_type==field_T) THEN 
    1213         nvert=6 
    1214          
    1215         DO ind=ind_b,ind_e 
    1216           d=>domain_type(ind) 
    1217          
    1218           DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    1219             DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    1220               IF (single .OR. d%assign_domain(i,j)==ind) ncell=ncell+1 
    1221             ENDDO 
    1222           ENDDO 
    1223  
    1224         END DO 
    1225        
    1226         status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
    1227         FieldId(NbField)=ncid 
     289     ENDIF 
     290      
     291     NbField=NbField+1 
     292     FieldName(NbField)=TRIM(ADJUSTL(name)) 
     293     FieldIndex(NbField)=1 
     294      
     295     PRINT *, 'Creating header for writefield : ', TRIM(FieldName(NbField)), ndomain, ind_b, ind_e ! FIXME 
     296     PRINT *, mesh_glo%ndomain, SIZE(mesh_glo%primal_own) 
     297     SELECT CASE(Field(ind_b)%field_type) 
     298     CASE(field_T) 
     299        PRINT *, '    field_type == field_T' ! FIXME 
     300        IF(single) THEN ! include halos 
     301           cells => mesh_glo%primal_all 
     302        ELSE 
     303           cells => mesh_glo%primal_own 
     304        END IF 
     305     CASE(field_Z) 
     306        PRINT *, '    field_type == field_Z' ! FIXME 
     307        IF(single) THEN ! include halos 
     308           cells => mesh_glo%dual_all 
     309        ELSE 
     310           cells => mesh_glo%dual_own 
     311        END IF 
     312     END SELECT 
     313      
     314     PRINT *, 'writefield : ind_b, ind_e :', ind_b, ind_e 
     315 
     316     ncell=0 
     317     DO ind=ind_b,ind_e 
     318        nvert=SIZE(cells(ind)%bnds_lon,1) 
     319        ncell = ncell + cells(ind)%ncell 
     320     END DO 
     321     PRINT *, 'writefield : found ',ncell, ' cells.' ! FIXME 
     322      
     323     status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
     324     FieldId(NbField)=ncid 
     325     status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
     326      
     327     level_size=0 
     328     SELECT CASE(Field(ind_b)%ndim) 
     329     CASE(2) 
     330        FieldVarId(NbField)%size=1 
     331        ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     332     CASE(3) 
     333        FieldVarId(NbField)%size=1 
     334        ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     335        status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id) 
     336        level_size=size(field(ind_b)%rval3d,2) 
     337     CASE(4) 
     338        FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
     339        ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
     340        status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id) 
     341        level_size=size(field(ind_b)%rval4d,2) 
     342        !          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id) 
     343     END SELECT 
     344      
     345     PRINT*,"create_header_gen : LEVEL_SIZE=",level_size 
     346     IF (level_size>0) THEN 
     347        status = NF90_DEF_VAR(ncid,'lev',NF90_DOUBLE,(/ dim3id /),levId) 
     348        status = NF90_PUT_ATT(ncid,levId,"axis","Z") 
     349     ENDIF 
     350      
     351     SELECT CASE(Field(ind_b)%field_type) 
     352     CASE(field_T) 
    1228353        status = NF90_DEF_DIM(ncid,'cell_i',ncell,ncellId) 
    1229354        status = NF90_DEF_DIM(ncid,'nvert_i',nvert,nvertid) 
    1230         level_size=0 
    1231         IF (Field(ind_b)%ndim==2)  THEN 
    1232           FieldVarId(NbField)%size=1 
    1233           ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    1234         ELSE IF (Field(ind_b)%ndim==3)  THEN 
    1235           FieldVarId(NbField)%size=1 
    1236           ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    1237           status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id) 
    1238           level_size=size(field(ind_b)%rval3d,2) 
    1239         ELSE IF (Field(1)%ndim==4) THEN 
    1240           FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
    1241           ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
    1242           status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id) 
    1243           level_size=size(field(ind_b)%rval4d,2) 
    1244 !          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id) 
    1245         ENDIF 
    1246         PRINT*,"LEVEL_SIZE=",level_size 
    1247  
    1248         status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
    1249         IF (level_size>0) THEN 
    1250           status = NF90_DEF_VAR(ncid,'lev',NF90_DOUBLE,(/ dim3id /),levId) 
    1251           status = NF90_PUT_ATT(ncid,levId,"axis","Z") 
    1252         ENDIF 
    1253          
    1254355        status = NF90_DEF_VAR(ncid,'lon_i',NF90_DOUBLE,(/ ncellId /),lonId) 
    1255         status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
    1256         status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
     356        status = NF90_DEF_VAR(ncid,'lat_i',NF90_DOUBLE,(/ ncellId /),latId) 
    1257357        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon_i") 
    1258         status = NF90_DEF_VAR(ncid,'lat_i',NF90_DOUBLE,(/ ncellId /),latId) 
    1259         status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
    1260         status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
    1261358        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat_i") 
    1262359        status = NF90_DEF_VAR(ncid,'bounds_lon_i',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
    1263360        status = NF90_DEF_VAR(ncid,'bounds_lat_i',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
    1264  
    1265         IF (Field(ind_b)%ndim==2) THEN 
    1266           IF (once) THEN  
    1267             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId /),FieldVarId(NbField)%nc_id(1)) 
    1268           ELSE 
    1269              status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
    1270           ENDIF 
    1271           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_i lat_i") 
    1272         ELSE IF (Field(ind_b)%ndim==3) THEN 
    1273           IF (once) THEN  
    1274             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(1)) 
    1275           ELSE 
    1276             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
    1277           ENDIF 
    1278           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_i lat_i") 
    1279         ELSE IF (Field(ind_b)%ndim==4) THEN 
    1280           DO i=1,FieldVarId(NbField)%size 
    1281             IF (once) THEN  
    1282               status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id /),  & 
    1283                                     FieldVarId(NbField)%nc_id(i)) 
    1284             ELSE 
    1285               status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /),  & 
    1286                                    FieldVarId(NbField)%nc_id(i)) 
    1287             ENDIF 
    1288             status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon_i lat_i") 
    1289           ENDDO         
    1290         ENDIF  
    1291    
    1292        
    1293         status = NF90_ENDDEF(ncid)       
    1294  
    1295         if (level_size>0) status = NF90_PUT_VAR(ncid,levId,(/ (l,l=1,level_size) /)) 
    1296           
    1297          ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell)) 
    1298            
    1299          n=0   
    1300          DO ind=ind_b,ind_e 
    1301            d=>domain_type(ind) 
    1302            DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    1303              DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    1304                IF (d%assign_domain(i,j)==ind .OR. single) THEN 
    1305                  n=n+1 
    1306                  CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n)) 
    1307                  lon(n)=lon(n)*180/Pi 
    1308                  lat(n)=lat(n)*180/Pi 
    1309                  DO k=0,5 
    1310                    CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n)) 
    1311                    bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
    1312                    bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
    1313                  ENDDO 
    1314                ENDIF 
    1315              ENDDO 
    1316            ENDDO 
    1317          ENDDO 
    1318          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /)) 
    1319          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /)) 
    1320          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
    1321          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
    1322    
    1323          DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
    1324  
    1325     ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
    1326         nvert=3 
    1327         DO ind=ind_b,ind_e 
    1328           d=>domain_type(ind) 
    1329          
    1330           DO j=d%jj_begin+1,d%jj_end 
    1331             DO i=d%ii_begin,d%ii_end-1 
    1332               ncell=ncell+1 
    1333             ENDDO 
    1334           ENDDO 
    1335  
    1336           DO j=d%jj_begin,d%jj_end-1 
    1337             DO i=d%ii_begin+1,d%ii_end 
    1338               ncell=ncell+1 
    1339             ENDDO 
    1340           ENDDO 
    1341  
    1342         END DO 
    1343        
    1344         status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
    1345         FieldId(NbField)=ncid 
     361     CASE(field_Z) 
    1346362        status = NF90_DEF_DIM(ncid,'cell_v',ncell,ncellId) 
    1347363        status = NF90_DEF_DIM(ncid,'nvert_v',nvert,nvertid) 
    1348  
    1349         IF (Field(ind_b)%ndim==2)  THEN 
    1350           FieldVarId(NbField)%size=1 
    1351           ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    1352         ELSE IF (Field(ind_b)%ndim==3)  THEN 
    1353           FieldVarId(NbField)%size=1 
    1354           ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    1355           status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id) 
    1356         ELSE IF (Field(1)%ndim==4) THEN 
    1357           FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
    1358           ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
    1359           status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id) 
    1360         ENDIF 
    1361  
    1362  
    1363        
    1364         status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
    1365        
    1366364        status = NF90_DEF_VAR(ncid,'lon_v',NF90_DOUBLE,(/ ncellId /),lonId) 
    1367         status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
    1368         status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
     365        status = NF90_DEF_VAR(ncid,'lat_v',NF90_DOUBLE,(/ ncellId /),latId) 
    1369366        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon_v") 
    1370         status = NF90_DEF_VAR(ncid,'lat_v',NF90_DOUBLE,(/ ncellId /),latId) 
    1371         status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
    1372         status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
    1373367        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat_v") 
    1374368        status = NF90_DEF_VAR(ncid,'bounds_lon_v',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
    1375369        status = NF90_DEF_VAR(ncid,'bounds_lat_v',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
    1376  
    1377  
    1378         IF (Field(ind_b)%ndim==2) THEN 
    1379           IF (once) THEN  
    1380             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId /),FieldVarId(NbField)%nc_id(1)) 
    1381           ELSE 
    1382             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
    1383           ENDIF 
    1384           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_v lat_v") 
    1385         ELSE IF (Field(ind_b)%ndim==3) THEN 
    1386           IF (once) THEN  
    1387             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(1)) 
    1388           ELSE 
    1389             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
    1390           ENDIF 
    1391           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_v lat_v") 
    1392         ELSE IF (Field(ind_b)%ndim==4) THEN 
    1393           DO q=1,FieldVarId(NbField)%size 
    1394             IF (once) THEN  
    1395               status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),ncprec,             & 
    1396                                     (/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(q)) 
    1397             ELSE 
    1398               status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),ncprec,             & 
    1399                                     (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q)) 
    1400             ENDIF 
    1401             status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon_v lat_v") 
    1402           ENDDO         
    1403         ENDIF  
    1404          
    1405         status = NF90_ENDDEF(ncid)       
    1406  
    1407          ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell)) 
    1408            
    1409          n=0   
    1410         
    1411          DO ind=ind_b,ind_e 
    1412            d=>domain_type(ind) 
    1413            DO j=d%jj_begin+1,d%jj_end 
    1414              DO i=d%ii_begin,d%ii_end-1 
    1415                nij=(j-1)*d%iim+i 
    1416                n=n+1 
    1417                CALL xyz2lonlat(d%vertex(:,vdown,i,j),lon(n),lat(n)) 
    1418                lon(n)=lon(n)*180/Pi 
    1419                lat(n)=lat(n)*180/Pi 
    1420  
    1421                CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n)) 
    1422                CALL xyz2lonlat(d%xyz(:,i,j-1),bounds_lon(1,n), bounds_lat(1,n)) 
    1423                CALL xyz2lonlat(d%xyz(:,i+1,j-1),bounds_lon(2,n), bounds_lat(2,n)) 
    1424  
    1425                DO k=0,2 
    1426                  bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
    1427                  bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
    1428                ENDDO 
    1429              ENDDO 
    1430            ENDDO 
    1431  
    1432            DO j=d%jj_begin,d%jj_end-1 
    1433              DO i=d%ii_begin+1,d%ii_end 
    1434                nij=(j-1)*d%iim+i 
    1435                n=n+1 
    1436                CALL xyz2lonlat(d%vertex(:,vup,i,j),lon(n),lat(n)) 
    1437                lon(n)=lon(n)*180/Pi 
    1438                lat(n)=lat(n)*180/Pi 
    1439                CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n)) 
    1440                CALL xyz2lonlat(d%xyz(:,i,j+1),bounds_lon(1,n), bounds_lat(1,n)) 
    1441                CALL xyz2lonlat(d%xyz(:,i-1,j+1),bounds_lon(2,n), bounds_lat(2,n)) 
    1442  
    1443                DO k=0,2 
    1444                  bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
    1445                  bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
    1446                ENDDO 
    1447              ENDDO 
    1448            ENDDO 
    1449          ENDDO  
    1450            
    1451          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /)) 
    1452          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /)) 
    1453          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
    1454          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
    1455    
    1456           DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
    1457       ENDIF 
    1458  
    1459  
    1460     END SUBROUTINE Create_Header_gen 
    1461  
    1462     SUBROUTINE Create_header_mpi(name,field,nind) 
    1463     USE netcdf_mod 
    1464     USE field_mod 
    1465     USE domain_mod 
    1466     USE spherical_geom_mod 
    1467     USE dimensions 
    1468     USE geometry 
    1469     USE mpi_mod 
    1470     USE mpipara 
    1471     IMPLICIT NONE 
    1472       CHARACTER(LEN=*) :: name 
    1473       CHARACTER(LEN=LEN_TRIM(ADJUSTL(name))) :: name_adj 
    1474       TYPE(t_field),POINTER :: field(:) 
    1475       INTEGER,OPTIONAL,INTENT(IN) :: nind 
    1476       INTEGER :: ncell 
    1477       INTEGER :: nvert 
    1478       REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:) 
    1479       TYPE(t_domain),POINTER :: d 
    1480       INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId 
    1481       INTEGER :: dim3id,dim4id 
    1482       INTEGER :: status 
    1483       INTEGER :: ind,i,j,k,n,q 
    1484       INTEGER :: iie,jje,iin,jjn 
    1485       INTEGER :: ind_b,ind_e 
    1486       INTEGER :: halo_size 
    1487       LOGICAL :: single  
    1488       INTEGER :: nij 
    1489       INTEGER :: ncell_glo(0:mpi_size-1) 
    1490       INTEGER :: displ, ncell_tot 
    1491        
    1492            
    1493       NbField=NbField+1 
    1494       name_adj=TRIM(ADJUSTL(name))  ! work around ICE with pgf90 
    1495       FieldName(NbField)=name_adj 
    1496       FieldIndex(NbField)=1 
    1497        
    1498       IF (PRESENT(nind)) THEN 
    1499         ind_b=nind 
    1500         ind_e=nind 
    1501         halo_size=1 
    1502         single=.TRUE. 
    1503       ELSE 
    1504         ind_b=1 
    1505         ind_e=ndomain 
    1506         halo_size=0 
    1507         single=.FALSE. 
    1508       ENDIF 
    1509        
    1510       ncell=0 
    1511        
    1512       IF (Field(ind_b)%field_type==field_T) THEN 
    1513         nvert=6 
    1514          
    1515         DO ind=ind_b,ind_e 
    1516           d=>domain(ind) 
    1517          
    1518           DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    1519             DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    1520               IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1 
    1521             ENDDO 
    1522           ENDDO 
    1523  
     370     END SELECT 
     371      
     372     status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
     373     status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
     374     status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
     375     status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
     376      
     377     SELECT CASE(Field(ind_b)%ndim) 
     378     CASE(2) 
     379        IF (once) THEN  
     380           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, & 
     381                (/ ncellId /),FieldVarId(NbField)%nc_id(1)) 
     382        ELSE 
     383           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, & 
     384                (/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
     385        END IF 
     386        CALL put_att_coordinates(1) 
     387     CASE(3) 
     388        IF (once) THEN  
     389           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, & 
     390                (/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(1)) 
     391        ELSE 
     392           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, & 
     393                (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
     394        END IF 
     395        CALL put_att_coordinates(1) 
     396     CASE(4) 
     397        DO i=1,FieldVarId(NbField)%size 
     398           IF (once) THEN  
     399              status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id /),  & 
     400                   FieldVarId(NbField)%nc_id(i)) 
     401           ELSE 
     402              status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /),  & 
     403                   FieldVarId(NbField)%nc_id(i)) 
     404           END IF 
     405           CALL put_att_coordinates(i) 
     406           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon_i lat_i") 
    1524407        END DO 
    1525        
    1526         CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr) 
    1527  
    1528         displ=0 
    1529         DO i=1,mpi_rank 
    1530           displ=displ+ncell_glo(i-1) 
    1531         ENDDO 
    1532         FieldVarId(NbField)%displ=displ 
    1533         ncell_tot=sum(ncell_glo(:)) 
    1534          
    1535 !        status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc', IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid) 
    1536         FieldId(NbField)=ncid 
    1537         status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId) 
    1538         status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
    1539  
    1540         IF (Field(ind_b)%ndim==2)  THEN 
    1541           FieldVarId(NbField)%size=1 
    1542           ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    1543         ELSE IF (Field(ind_b)%ndim==3)  THEN 
    1544           FieldVarId(NbField)%size=1 
    1545           ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    1546           status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id) 
    1547         ELSE IF (Field(1)%ndim==4) THEN 
    1548           FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
    1549           ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
    1550           status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id) 
    1551 !          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id) 
    1552         ENDIF 
    1553        
    1554         status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
    1555        
    1556         status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) 
    1557         status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
    1558         status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
    1559         status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") 
    1560         status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) 
    1561         status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
    1562         status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
    1563         status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") 
    1564         status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
    1565         status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
    1566  
    1567         IF (Field(ind_b)%ndim==2) THEN 
    1568           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
    1569           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
    1570           status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/)) 
    1571         ELSE IF (Field(ind_b)%ndim==3) THEN 
    1572           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
    1573           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
    1574           status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED,   & 
    1575                                          (/ncell_tot,size(field(ind_b)%rval3d,2),1/)) 
    1576         ELSE IF (Field(ind_b)%ndim==4) THEN 
    1577           DO i=1,FieldVarId(NbField)%size 
    1578             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /),  & 
    1579                                   FieldVarId(NbField)%nc_id(i)) 
    1580             status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat") 
    1581             status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED,   & 
    1582                                            (/ncell_tot,size(field(ind_b)%rval4d,2),1/)) 
    1583           ENDDO         
    1584         ENDIF  
    1585    
    1586        
    1587         status = NF90_ENDDEF(ncid)       
    1588  
    1589         ncell=1 
    1590         DO ind=ind_b,ind_e 
    1591           d=>domain(ind) 
    1592    
    1593           n=0 
    1594           DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    1595             DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    1596               IF (single .OR. d%own(i,j)) n=n+1 
    1597             ENDDO 
    1598           ENDDO 
    1599            
    1600          ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n)) 
    1601            
    1602           n=0   
    1603           DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    1604             DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    1605                 IF (d%own(i,j) .OR. single) THEN 
    1606                 n=n+1 
    1607                 CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n)) 
    1608                 lon(n)=lon(n)*180/Pi 
    1609                 lat(n)=lat(n)*180/Pi 
    1610                 DO k=0,5 
    1611                   CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n)) 
    1612                   bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
    1613                   bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
    1614                 ENDDO 
    1615               ENDIF 
    1616             ENDDO 
    1617           ENDDO 
    1618           status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /)) 
    1619           status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /)) 
    1620           status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) 
    1621           status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) 
    1622    
    1623           ncell=ncell+n 
    1624           DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
    1625       END DO  
    1626  
    1627     ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
    1628         nvert=3 
    1629         DO ind=ind_b,ind_e 
    1630           d=>domain(ind) 
    1631          
    1632           DO j=d%jj_begin+1,d%jj_end 
    1633             DO i=d%ii_begin,d%ii_end-1 
    1634               ncell=ncell+1 
    1635             ENDDO 
    1636           ENDDO 
    1637  
    1638           DO j=d%jj_begin,d%jj_end-1 
    1639             DO i=d%ii_begin+1,d%ii_end 
    1640               ncell=ncell+1 
    1641             ENDDO 
    1642           ENDDO 
    1643  
    1644         END DO 
    1645          
    1646         CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr) 
    1647  
    1648         displ=0 
    1649         DO i=1,mpi_rank 
    1650           displ=displ+ncell_glo(i-1) 
    1651         ENDDO 
    1652         FieldVarId(NbField)%displ=displ 
    1653         ncell_tot=sum(ncell_glo(:)) 
    1654                
    1655 !        status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc',IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid) 
    1656 !        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
    1657         FieldId(NbField)=ncid 
    1658         status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId) 
    1659         status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
    1660  
    1661         IF (Field(ind_b)%ndim==2)  THEN 
    1662           FieldVarId(NbField)%size=1 
    1663           ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    1664         ELSE IF (Field(ind_b)%ndim==3)  THEN 
    1665           FieldVarId(NbField)%size=1 
    1666           ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    1667           status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id) 
    1668         ELSE IF (Field(1)%ndim==4) THEN 
    1669           FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
    1670           ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
    1671           status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id) 
    1672         ENDIF 
    1673  
    1674  
    1675        
    1676         status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
    1677        
    1678         status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) 
    1679         status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
    1680         status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
    1681         status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") 
    1682         status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) 
    1683         status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
    1684         status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
    1685         status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") 
    1686         status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
    1687         status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
    1688  
    1689  
    1690         IF (Field(ind_b)%ndim==2) THEN 
    1691           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
    1692           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
    1693           status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/)) 
    1694         ELSE IF (Field(ind_b)%ndim==3) THEN 
    1695           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
    1696           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
    1697           status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED,   & 
    1698                                          (/ncell_tot,size(field(ind_b)%rval3d,2),1/)) 
    1699         ELSE IF (Field(ind_b)%ndim==4) THEN 
    1700           DO q=1,FieldVarId(NbField)%size 
    1701             status = NF90_DEF_VAR(ncid,name_adj//int2str(q),ncprec,             & 
    1702                                   (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q)) 
    1703             status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat") 
    1704            status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED, & 
    1705                                           (/ncell_tot,size(field(ind_b)%rval4d,2),1/)) 
    1706           ENDDO         
    1707         ENDIF  
    1708          
    1709         status = NF90_ENDDEF(ncid)       
    1710  
    1711         ncell=1 
    1712         DO ind=ind_b,ind_e 
    1713           d=>domain(ind) 
    1714           CALL swap_geometry(ind) 
    1715           CALL swap_dimensions(ind) 
    1716    
    1717           n=0 
    1718           DO j=jj_begin+1,jj_end 
    1719             DO i=ii_begin,ii_end-1 
    1720               n=n+1 
    1721             ENDDO 
    1722           ENDDO 
    1723  
    1724           DO j=jj_begin,jj_end-1 
    1725             DO i=ii_begin+1,ii_end 
    1726               n=n+1 
    1727             ENDDO 
    1728           ENDDO 
    1729  
    1730          ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n)) 
    1731            
    1732           n=0   
    1733         
    1734           DO j=jj_begin+1,jj_end 
    1735             DO i=ii_begin,ii_end-1 
    1736               nij=(j-1)*iim+i 
    1737               n=n+1 
    1738               CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n)) 
    1739               lon(n)=lon(n)*180/Pi 
    1740               lat(n)=lat(n)*180/Pi 
    1741               CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
    1742               CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) 
    1743               CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) 
    1744  
    1745               DO k=0,2 
    1746                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
    1747                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
    1748               ENDDO 
    1749             ENDDO 
    1750           ENDDO 
    1751  
    1752           DO j=jj_begin,jj_end-1 
    1753             DO i=ii_begin+1,ii_end 
    1754               nij=(j-1)*iim+i 
    1755               n=n+1 
    1756               CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n)) 
    1757               lon(n)=lon(n)*180/Pi 
    1758               lat(n)=lat(n)*180/Pi 
    1759               CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
    1760               CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) 
    1761               CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) 
    1762  
    1763               DO k=0,2 
    1764                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
    1765                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
    1766               ENDDO 
    1767             ENDDO 
    1768           ENDDO 
    1769            
    1770            
    1771           status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /)) 
    1772           status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /)) 
    1773           status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) 
    1774           status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) 
    1775    
    1776           ncell=ncell+n 
    1777           DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
    1778       END DO           
    1779     ENDIF 
    1780  
    1781  
    1782    END SUBROUTINE Create_Header_mpi   
    1783     
     408     END SELECT 
     409      
     410     status = NF90_ENDDEF(ncid)       
     411      
     412     if (level_size>0) status = NF90_PUT_VAR(ncid,levId,(/ (l,l=1,level_size) /)) 
     413      
     414     ALLOCATE(lon(ncell),lat(ncell)) 
     415     ALLOCATE(bounds_lon(nvert,ncell), bounds_lat(nvert,ncell)) 
     416     n_begin=0   
     417     DO ind=ind_b,ind_e 
     418        n_end = n_begin + cells(ind)%ncell 
     419        PRINT *, 'create_header_gen ', n_begin, n_end, SHAPE(cells(ind)%bnds_lon), SHAPE(bounds_lon) 
     420        lat(n_begin+1:n_end) = cells(ind)%lat(:) *180./Pi 
     421        lon(n_begin+1:n_end) = cells(ind)%lon(:) *180./Pi 
     422        bounds_lon(:,n_begin+1:n_end) = cells(ind)%bnds_lon(:,:) *180./Pi 
     423        bounds_lat(:,n_begin+1:n_end) = cells(ind)%bnds_lat(:,:) *180./Pi 
     424        n_begin = n_end 
     425     END DO 
     426     status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /)) 
     427     status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /)) 
     428     status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
     429     status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
     430      
     431     DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
     432      
     433   CONTAINS 
     434      
     435     SUBROUTINE put_att_coordinates(ind) 
     436       INTEGER :: ind 
     437       SELECT CASE(Field(ind_b)%field_type) 
     438       CASE(field_T) 
     439          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(ind), & 
     440               "coordinates","lon_i lat_i") 
     441       CASE(field_Z) 
     442          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(ind), & 
     443               "coordinates","lon_v lat_v")       
     444       END SELECT 
     445     END SUBROUTINE put_att_coordinates 
     446      
     447   END SUBROUTINE Create_Header_gen 
     448     
    1784449   SUBROUTINE Close_files 
    1785450   USE netcdf_mod 
    1786    IMPLICIT NONE 
    1787451     INTEGER :: i,k,status 
    1788452!$OMP MASTER      
    … …  
    1794458      
    1795459      
    1796   function int2str(int) 
    1797     implicit none 
    1798     integer, parameter :: MaxLen=10 
    1799     integer,intent(in) :: int 
    1800     character(len=MaxLen) :: int2str 
    1801     logical :: flag 
    1802     integer :: i 
    1803     flag=.true. 
    1804      
    1805     i=int 
    1806      
    1807     int2str='' 
    1808     do while (flag) 
    1809       int2str=CHAR(MOD(i,10)+48)//int2str 
    1810       i=i/10 
    1811       if (i==0) flag=.false. 
    1812     enddo 
    1813   end function int2str 
    1814  
    1815460end module write_field_mod 
    1816    
Note: See TracChangeset for help on using the changeset viewer.