Changeset 880


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

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

Location:
codes/icosagcm/devel/src
Files:
3 added
3 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    

  • codes/icosagcm/devel/src/output/xios_mod.F90

    r874 r880  
    11MODULE xios_mod 
    2  
     2   
    33#ifdef CPP_USING_XIOS 
    4    USE xios 
     4  USE xios 
    55#endif 
    6    IMPLICIT NONE  
    7  
    8   PUBLIC 
    9   LOGICAL,SAVE :: using_xios 
    10  
    11   INTEGER,SAVE :: ncell_i 
    12 !$OMP THREADPRIVATE(ncell_i) 
    13   INTEGER,SAVE :: ncell_v 
    14 !$OMP THREADPRIVATE(ncell_v) 
    15   INTEGER,SAVE :: ncell_e 
    16 !$OMP THREADPRIVATE(ncell_e) 
    17  
    18   PRIVATE ncell_i,ncell_v,ncell_e 
     6 
     7  USE prec, ONLY : rstd 
     8  USE field_mod, ONLY : t_field, field_T, field_U, field_Z 
     9  USE domain_mod, ONLY : t_domain, t_cellset, domain, ndomain, mesh_loc 
     10  IMPLICIT NONE    
     11  PRIVATE 
     12  SAVE 
     13 
     14  LOGICAL :: using_xios 
     15 
     16  INTEGER :: ncell_i, ncell_v, ncell_e 
     17!$OMP THREADPRIVATE(ncell_i, ncell_v, ncell_e) 
     18 
     19  PUBLIC :: using_xios, xios_init, & 
     20       xios_init_write_field,  xios_write_field_finalize, & 
     21       xios_write_field, xios_read_field 
    1922 
    2023#ifdef CPP_USING_XIOS 
     24 
     25  PUBLIC :: xios_timestep,    &  
     26       xios_set_file_attr, xios_set_fieldgroup_attr, & 
     27       xios_set_filegroup_attr, xios_get_axis_attr, & 
     28       xios_send_field, xios_read_var, & 
     29       xios_update_calendar, xios_set_context 
    2130   
    2231CONTAINS 
    2332  
    2433 SUBROUTINE xios_init 
    25    USE getin_mod 
    26    USE xios 
    27    USE mpipara 
    28    IMPLICIT NONE 
    29     TYPE(xios_context) :: ctx_hdl 
    30  
    31      using_xios=.TRUE. 
    32      CALL xios_context_initialize("icosagcm",comm_icosa) 
    33      CALL xios_get_handle("icosagcm",ctx_hdl) 
    34      CALL xios_set_current_context(ctx_hdl) 
    35     
    36  END SUBROUTINE xios_init   
    37   
     34   USE mpipara, ONLY : comm_icosa 
     35   TYPE(xios_context) :: ctx_hdl 
     36 
     37   using_xios=.TRUE. 
     38   CALL xios_context_initialize("icosagcm",comm_icosa) 
     39   CALL xios_get_handle("icosagcm",ctx_hdl) 
     40   CALL xios_set_current_context(ctx_hdl)    
     41 
     42 END SUBROUTINE xios_init 
     43 
    3844 SUBROUTINE xios_init_write_field 
    39  USE genmod 
    40  USE mpipara 
    41  USE xios 
    42  USE grid_param 
    43  USE domain_mod 
    44  USE dimensions 
    45  USE spherical_geom_mod 
    46  USE geometry 
    47  USE mpi_mod 
    48  USE time_mod 
    49  USE metric, ONLY : vup,vdown, cell_glo 
    50  USE disvert_mod, ONLY : presnivs 
    51  IMPLICIT NONE 
    52   TYPE(xios_context) :: ctx_hdl 
    53   TYPE(xios_duration)      :: dtime 
    54   REAL(rstd) :: lev_value(llm) 
    55   REAL(rstd) :: lev_valuep1(llm+1) 
    56   REAL(rstd) :: nq_value(nqtot) 
    57   INTEGER :: ncell, ncell_tot, ncell_glo(0:mpi_size-1), displ 
    58   INTEGER :: ind, i,j,k,l,ij 
    59   REAL(rstd),ALLOCATABLE    :: lon(:), lat(:), bounds_lon(:,:), bounds_lat(:,:) 
    60   INTEGER, ALLOCATABLE      :: ind_glo(:) 
    61   TYPE(t_domain),POINTER :: d 
    62  
    63 !$OMP BARRIER 
    64 !$OMP MASTER 
    65 !   CALL xios_context_initialize("icosagcm",comm_icosa) 
     45   USE time_mod, ONLY : dt, itau_out 
     46   USE grid_param, ONLY : llm, nqtot 
     47   USE mpi_mod, ONLY : MPI_INTEGER 
     48 
     49   TYPE(xios_context) :: ctx_hdl 
     50   TYPE(xios_duration)      :: dtime 
     51   REAL(rstd) :: lev_value(llm) 
     52   REAL(rstd) :: lev_valuep1(llm+1) 
     53   REAL(rstd) :: nq_value(nqtot) 
     54   INTEGER :: l, ncell, ncell_tot, displ 
     55   REAL(rstd),ALLOCATABLE    :: lon(:), lat(:), bounds_lon(:,:), bounds_lat(:,:) 
     56   INTEGER, ALLOCATABLE      :: ind_glo(:) 
     57   TYPE(t_domain),POINTER :: d 
     58    
     59   !$OMP BARRIER 
     60   !$OMP MASTER 
     61   !   CALL xios_context_initialize("icosagcm",comm_icosa) 
    6662   CALL xios_get_handle("icosagcm",ctx_hdl) 
    6763   CALL xios_set_current_context(ctx_hdl) 
    … …  
    7167   CALL xios_set_axis_attr("lev",n_glo=llm ,value=lev_value) ; 
    7268   CALL xios_set_axis_attr("levp1",n_glo=llm+1 ,value=lev_valuep1) ; 
    73    CALL xios_set_axis_attr("presnivs_mb",n_glo=llm, value=presnivs/100., unit="mb") ;  
    7469   CALL xios_set_axis_attr("nq",n_glo=nqtot, value=nq_value) ; 
    7570    
    76    ncell=0 
    77    DO ind=1,ndomain 
    78      d=>domain(ind) 
    79          
    80      DO j=d%jj_begin,d%jj_end 
    81        DO i=d%ii_begin,d%ii_end 
    82          IF (domain(ind)%own(i,j)) ncell=ncell+1 
    83        ENDDO 
    84      ENDDO 
    85    ENDDO       
    86    ncell_i=ncell 
    87     
    88    CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr) 
    89  
    90    displ=0 
    91    DO i=1,mpi_rank 
    92      displ=displ+ncell_glo(i-1) 
    93    ENDDO 
    94  
    95    ncell_tot=sum(ncell_glo(:)) 
    96     
    97    ALLOCATE(lon(ncell), lat(ncell), bounds_lon(0:5,ncell), bounds_lat(0:5,ncell), ind_glo(ncell))  
    98     
    99    ncell=0 
    100    DO ind=1,ndomain 
    101      d=>domain(ind) 
    102          
    103      DO j=d%jj_begin,d%jj_end 
    104        DO i=d%ii_begin,d%ii_end 
    105          IF (domain(ind)%own(i,j)) THEN  
    106            ncell=ncell+1 
    107            CALL xyz2lonlat(d%xyz(:,i,j),lon(ncell),lat(ncell)) 
    108            lon(ncell)=lon(ncell)*180/Pi 
    109            lat(ncell)=lat(ncell)*180/Pi 
    110            DO k=0,5 
    111              CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,ncell), bounds_lat(k,ncell)) 
    112              bounds_lat(k,ncell)=bounds_lat(k,ncell)*180/Pi 
    113              bounds_lon(k,ncell)=bounds_lon(k,ncell)*180/Pi 
    114            ENDDO 
    115            ind_glo(ncell)=domain(ind)%assign_cell_glo(i,j)-1  
    116          ENDIF 
    117        ENDDO 
    118      ENDDO 
    119    ENDDO          
    120  
     71   !------------------ primal cells ------------------ 
     72   CALL collect_bounds(6, mesh_loc%primal_own) 
     73   ncell_i=ncell    
    12174   CALL xios_set_domaingroup_attr("i",ni_glo=ncell_tot, ibegin=displ, ni=ncell) 
    12275   CALL xios_set_domaingroup_attr("i", data_dim=1, type='unstructured' , nvertex=6, i_index=ind_glo) 
    12376   CALL xios_set_domaingroup_attr("i",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat) 
    124     
    12577   DEALLOCATE(lon, lat, bounds_lon, bounds_lat,ind_glo)  
    126      
    127  
    128  
    129    ncell=0 
    130    DO ind=1,ndomain 
    131      d=>domain(ind) 
    132  
    133      DO j=d%jj_begin,d%jj_end 
    134        DO i=d%ii_begin,d%ii_end 
    135          DO k=0,5 
    136             IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j & 
    137                                                                                  .AND. d%edge_assign_pos(k,i,j)==k) THEN 
    138                ncell=ncell+1 
    139             ENDIF 
    140          ENDDO 
    141        ENDDO 
    142      ENDDO 
    143    ENDDO 
     78    
     79   !--------------------- dual cells ------------------ 
     80    
     81   CALL collect_bounds(3, mesh_loc%dual_own) 
     82   ncell_v=ncell 
     83   CALL xios_set_domain_attr("v",ni_glo=ncell_tot, ibegin=displ, ni=ncell) 
     84   CALL xios_set_domain_attr("v", data_dim=1, type='unstructured' , nvertex=3) 
     85   CALL xios_set_domain_attr("v",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat) 
     86   DEALLOCATE(lon, lat, bounds_lon, bounds_lat,ind_glo)  
     87    
     88   !---------------------- edges ----------------------- 
     89    
     90   CALL collect_bounds(2, mesh_loc%edge_own) 
    14491   ncell_e=ncell 
    145     
    146    CALL MPI_ALLGATHER(ncell_e,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr) 
    147    displ=0 
    148    DO i=1,mpi_rank 
    149      displ=displ+ncell_glo(i-1) 
    150    ENDDO 
    151    ncell_tot=sum(ncell_glo(:)) 
    152     
    153    ALLOCATE(lon(ncell), lat(ncell), bounds_lon(0:1,ncell), bounds_lat(0:1,ncell),ind_glo(ncell))  
    154  
    155  
    156    ncell=0 
    157    DO ind=1,ndomain 
    158      d=>domain(ind) 
    159      CALL swap_dimensions(ind) 
    160      CALL swap_geometry(ind) 
    161  
    162      DO j=d%jj_begin,d%jj_end 
    163        DO i=d%ii_begin,d%ii_end 
    164          DO k=0,5 
    165            IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j & 
    166                                                                                  .AND. d%edge_assign_pos(k,i,j)==k) THEN 
    167               ncell=ncell+1 
    168               ij=(j-1)*iim+i 
    169  
    170               lon(ncell)=lon_e(ij+u_pos(k+1))*180/Pi 
    171               lat(ncell)=lat_e(ij+u_pos(k+1))*180/Pi 
    172                 
    173               CALL xyz2lonlat(d%vertex(:,MOD((k-1)+6,6),i,j),bounds_lon(0,ncell), bounds_lat(0,ncell)) 
    174               CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(1,ncell), bounds_lat(1,ncell)) 
    175               bounds_lon(:,ncell)=bounds_lon(:,ncell)*180/Pi 
    176               bounds_lat(:,ncell)=bounds_lat(:,ncell)*180/Pi 
    177               ind_glo(ncell)=cell_glo(d%assign_cell_glo(i,j))%edge(MOD(k+d%delta(i,j)+6,6))-1  
    178            ENDIF                
    179          ENDDO 
    180        ENDDO 
    181      ENDDO 
    182    ENDDO 
    18392   CALL xios_set_domain_attr("u",ni_glo=ncell_tot, ibegin=displ, ni=ncell) 
    18493   CALL xios_set_domain_attr("u", data_dim=1, type='unstructured' , nvertex=2, i_index=ind_glo) 
    18594   CALL xios_set_domain_attr("u",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat) 
    186  
    18795   DEALLOCATE(lon, lat, bounds_lon, bounds_lat,ind_glo)  
    188  
    189  
    190    ncell=0 
    191    DO ind=1,ndomain 
    192      d=>domain(ind) 
    193          
    194      DO j=d%jj_begin+1,d%jj_end 
    195        DO i=d%ii_begin,d%ii_end-1 
    196          ncell=ncell+1 
    197        ENDDO 
    198      ENDDO 
    199  
    200      DO j=d%jj_begin,d%jj_end-1 
    201        DO i=d%ii_begin+1,d%ii_end 
    202           ncell=ncell+1 
    203         ENDDO 
    204      ENDDO 
    205  
    206    ENDDO       
    207     
    208    ncell_v=ncell 
    209     
    210    CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr) 
    211  
    212    displ=0 
    213    DO i=1,mpi_rank 
    214      displ=displ+ncell_glo(i-1) 
    215    ENDDO 
    216  
    217    ncell_tot=sum(ncell_glo(:)) 
    218     
    219    ALLOCATE(lon(ncell), lat(ncell), bounds_lon(0:2,ncell), bounds_lat(0:2,ncell))  
    220     
    221    ncell=0 
    222    DO ind=1,ndomain 
    223      d=>domain(ind) 
    224   
    225      DO j=d%jj_begin+1,d%jj_end 
    226        DO i=d%ii_begin,d%ii_end-1 
    227            ncell=ncell+1 
    228            CALL xyz2lonlat(d%vertex(:,vdown,i,j),lon(ncell),lat(ncell)) 
    229            lon(ncell)=lon(ncell)*180/Pi 
    230            lat(ncell)=lat(ncell)*180/Pi 
    231  
    232            CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,ncell), bounds_lat(0,ncell)) 
    233            CALL xyz2lonlat(d%xyz(:,i,j-1),bounds_lon(1,ncell), bounds_lat(1,ncell)) 
    234            CALL xyz2lonlat(d%xyz(:,i+1,j-1),bounds_lon(2,ncell), bounds_lat(2,ncell)) 
    235  
    236            DO k=0,2 
    237              bounds_lat(k,ncell)=bounds_lat(k,ncell)*180/Pi 
    238              bounds_lon(k,ncell)=bounds_lon(k,ncell)*180/Pi 
    239            ENDDO 
    240          ENDDO 
    241        ENDDO  
    242   
    243        DO j=d%jj_begin,d%jj_end-1 
    244          DO i=d%ii_begin+1,d%ii_end 
    245            ncell=ncell+1 
    246            CALL xyz2lonlat(d%vertex(:,vup,i,j),lon(ncell),lat(ncell)) 
    247            lon(ncell)=lon(ncell)*180/Pi 
    248            lat(ncell)=lat(ncell)*180/Pi 
    249            CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,ncell), bounds_lat(0,ncell)) 
    250            CALL xyz2lonlat(d%xyz(:,i,j+1),bounds_lon(1,ncell), bounds_lat(1,ncell)) 
    251            CALL xyz2lonlat(d%xyz(:,i-1,j+1),bounds_lon(2,ncell), bounds_lat(2,ncell)) 
    252  
    253            DO k=0,2 
    254              bounds_lat(k,ncell)=bounds_lat(k,ncell)*180/Pi 
    255              bounds_lon(k,ncell)=bounds_lon(k,ncell)*180/Pi 
    256            ENDDO 
    257          ENDDO 
    258        ENDDO  
    259          
    260    ENDDO          
    261  
    262    
    263    CALL xios_set_domain_attr("v",ni_glo=ncell_tot, ibegin=displ, ni=ncell) 
    264    CALL xios_set_domain_attr("v", data_dim=1, type='unstructured' , nvertex=3) 
    265    CALL xios_set_domain_attr("v",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat) 
    266  
    267  
     96    
    26897   dtime%second=dt 
    26998   CALL xios_set_timestep(dtime) 
    27099    
    271100   CALL xios_set_fieldgroup_attr("standard_output", freq_op=itau_out*xios_timestep, freq_offset=(itau_out-1)*xios_timestep) 
    272      
     101    
    273102   CALL xios_close_context_definition() 
    274 !$OMP END MASTER 
    275 !$OMP BARRIER 
     103   !$OMP END MASTER 
     104   !$OMP BARRIER 
     105    
     106 CONTAINS 
     107    
     108   SUBROUTINE collect_bounds(nvert, cells) 
     109     USE mpipara, ONLY : comm_icosa, mpi_size, mpi_rank 
     110     INTEGER, INTENT(IN) :: nvert 
     111     TYPE(t_cellset)     :: cells(:) 
     112     INTEGER :: i, ind, n_beg, n_end, ierr, ncell_glo(0:mpi_size-1) 
     113     ncell = SUM(cells%ncell) 
     114     CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER, & 
     115          ncell_glo,1, MPI_INTEGER, comm_icosa, ierr) 
     116     displ=0 
     117     DO i=1,mpi_rank 
     118        displ=displ+ncell_glo(i-1) 
     119     ENDDO 
     120     ncell_tot=sum(ncell_glo(:)) 
     121      
     122     ALLOCATE(lon(ncell), lat(ncell), ind_glo(ncell)) 
     123     ALLOCATE(bounds_lon(nvert,ncell), bounds_lat(nvert,ncell)) 
     124      
     125     n_beg=0 
     126     DO ind=1,ndomain 
     127        n_end = n_beg + cells(ind)%ncell 
     128        ind_glo(n_beg+1:n_end) = cells(ind)%ind_glo(:) 
     129        lon(n_beg+1:n_end) = cells(ind)%lon(:) 
     130        lat(n_beg+1:n_end) = cells(ind)%lat(:) 
     131        bounds_lon(:,n_beg+1:n_end) = cells(ind)%bnds_lon(:,:) 
     132        bounds_lat(:,n_beg+1:n_end) = cells(ind)%bnds_lat(:,:) 
     133        n_beg = n_end 
     134     END DO 
     135   END SUBROUTINE collect_bounds 
    276136    
    277137 END SUBROUTINE xios_init_write_field 
    278    
    279    
     138  
    280139 SUBROUTINE xios_write_field(name,field) 
    281  USE field_mod 
    282  IMPLICIT NONE  
    283140   CHARACTER(LEN=*),INTENT(IN) :: name 
    284141   TYPE(t_field), POINTER :: field(:) 
    285    CHARACTER(LEN=10) :: str_number 
    286    INTEGER :: iq 
     142   TYPE(t_cellset), POINTER :: cells(:) 
     143   INTEGER :: ncells 
    287144 
    288145!$OMP BARRIER 
    289146!$OMP MASTER 
    290     
    291    IF (Field(1)%field_type==field_T) THEN 
    292      IF (field(1)%ndim==2) THEN 
    293         CALL xios_write_field_scalar(name,field,1,1) 
    294      ELSE IF (field(1)%ndim==3) THEN 
    295         CALL xios_write_field_scalar(name,field,size(field(1)%rval3d,2),1) 
    296      ELSE IF (field(1)%ndim==4) THEN 
    297 !        DO iq=1,size(field(1)%rval4d,3) 
    298 !          WRITE(str_number,'(i10)') iq 
    299 !          CALL xios_write_field_scalar(name//TRIM(ADJUSTL(str_number)),field,size(field(1)%rval4d,2),iq) 
    300           CALL xios_write_field_scalar(name,field,size(field(1)%rval4d,2),size(field(1)%rval4d,3)) 
    301 !        ENDDO 
    302      ELSE 
    303         PRINT *, "xios_write_field : dimension > 4 are not supported for now" 
    304      ENDIF 
    305  
    306    ELSE IF (Field(1)%field_type==field_U) THEN 
    307       IF (field(1)%ndim==2) THEN 
    308         CALL xios_write_field_U(name,field,1,1) 
    309       ELSE IF (field(1)%ndim==3) THEN 
    310         CALL xios_write_field_U(name,field,size(field(1)%rval3d,2),1) 
    311       ELSE IF (field(1)%ndim==4) THEN 
    312         CALL xios_write_field_U(name,field,size(field(1)%rval4d,2),size(field(1)%rval4d,3)) 
    313       ELSE 
    314         PRINT *, "xios_write_field : dimension > 4 are not supported for now" 
    315       ENDIF 
    316  
    317     ELSE IF (Field(1)%field_type==field_Z) THEN 
    318      IF (field(1)%ndim==2) THEN 
    319         CALL xios_write_field_vort(name,field,1) 
    320       ELSE IF (field(1)%ndim==3) THEN 
    321         CALL xios_write_field_vort(name,field,size(field(1)%rval3d,2)) 
    322       ELSE IF (field(1)%ndim==4) THEN 
    323         DO iq=1,size(field(1)%rval4d,3) 
    324           WRITE(str_number,'(i10)') iq 
    325           CALL xios_write_field_vort(name//TRIM(ADJUSTL(str_number)),field,size(field(1)%rval4d,2),iq) 
    326         ENDDO 
    327       ELSE 
    328         PRINT *, "xios_write_field : dimension > 4 are not supported for now" 
    329       ENDIF 
    330     ENDIF 
     147 
     148   SELECT CASE(field(1)%field_type) 
     149   CASE(field_T) 
     150      cells => mesh_loc%primal_own 
     151      ncells = ncell_i 
     152   CASE(field_U) 
     153      cells => mesh_loc%edge_own 
     154      ncells = ncell_e 
     155   CASE(field_Z) 
     156      cells => mesh_loc%dual_own 
     157      ncells = ncell_v 
     158   END SELECT 
     159 
     160   IF (field(1)%ndim>4) THEN 
     161      PRINT *, "xios_write_field : dimension > 4 are not supported for now" 
     162   ELSE 
     163      CALL xios_write_field_hex(name, field, cells, &  
     164           ncells, field(1)%dim3, field(1)%dim4) 
     165   END IF 
     166    
    331167!$OMP END MASTER 
    332168!$OMP BARRIER 
    … …  
    335171 
    336172 SUBROUTINE xios_read_field(name,field) 
    337  USE field_mod 
    338  IMPLICIT NONE  
    339173   CHARACTER(LEN=*),INTENT(IN) :: name 
    340174   TYPE(t_field), POINTER :: field(:) 
    341    CHARACTER(LEN=10) :: str_number 
    342    INTEGER :: iq 
     175   TYPE(t_cellset), POINTER :: cells(:) 
     176   INTEGER :: ncells 
    343177 
    344178!$OMP BARRIER 
    345179!$OMP MASTER 
    346     
    347    IF (Field(1)%field_type==field_T) THEN 
    348      IF (field(1)%ndim==2) THEN 
    349         CALL xios_read_field_scalar(name,field,1,1) 
    350       ELSE IF (field(1)%ndim==3) THEN 
    351         CALL xios_read_field_scalar(name,field,size(field(1)%rval3d,2),1) 
    352       ELSE IF (field(1)%ndim==4) THEN 
    353 !        DO iq=1,size(field(1)%rval4d,3) 
    354 !          WRITE(str_number,'(i10)') iq 
    355 !          CALL xios_read_field_scalar(name//TRIM(ADJUSTL(str_number)),field,size(field(1)%rval4d,2),iq) 
    356 !        ENDDO 
    357           CALL xios_read_field_scalar(name,field,size(field(1)%rval4d,2),size(field(1)%rval4d,3)) 
    358       ELSE 
    359         PRINT *, "xios_write_field : dimension > 4 are not supported for now" 
    360       ENDIF 
    361    ELSE IF (Field(1)%field_type==field_U) THEN 
    362      IF (field(1)%ndim==2) THEN 
    363         CALL xios_read_field_u(name,field,1,1) 
    364       ELSE IF (field(1)%ndim==3) THEN 
    365         CALL xios_read_field_u(name,field,size(field(1)%rval3d,2),1) 
    366       ELSE IF (field(1)%ndim==4) THEN 
    367           CALL xios_read_field_u(name,field,size(field(1)%rval4d,2),size(field(1)%rval4d,3)) 
    368       ELSE 
    369         PRINT *, "xios_write_field : dimension > 4 are not supported for now" 
    370       ENDIF 
    371     ELSE IF (Field(1)%field_type==field_Z) THEN 
    372      IF (field(1)%ndim==2) THEN 
    373         CALL xios_read_field_vort(name,field,1) 
    374       ELSE IF (field(1)%ndim==3) THEN 
    375         CALL xios_read_field_vort(name,field,size(field(1)%rval3d,2)) 
    376       ELSE IF (field(1)%ndim==4) THEN 
    377         DO iq=1,size(field(1)%rval4d,3) 
    378           WRITE(str_number,'(i10)') iq 
    379           CALL xios_read_field_vort(name//TRIM(ADJUSTL(str_number)),field,size(field(1)%rval4d,2),iq) 
    380         ENDDO 
    381       ELSE 
    382         PRINT *, "xios_write_field : dimension > 4 are not supported for now" 
    383       ENDIF 
    384     ENDIF 
     180 
     181   SELECT CASE(field(1)%field_type) 
     182   CASE(field_T) 
     183      cells => mesh_loc%primal_own 
     184      ncells = ncell_i 
     185   CASE(field_U) 
     186      cells => mesh_loc%edge_own 
     187      ncells = ncell_e 
     188   CASE(field_Z) 
     189      cells => mesh_loc%dual_own 
     190      ncells = ncell_v 
     191   END SELECT 
     192 
     193   IF (field(1)%ndim>4) THEN 
     194      PRINT *, "xios_read_field : dimension > 4 are not supported for now" 
     195   ELSE 
     196      CALL xios_read_field_hex(name, field, cells, &  
     197           ncells, field(1)%dim3, field(1)%dim4) 
     198   END IF 
     199    
    385200!$OMP END MASTER 
    386201!$OMP BARRIER 
    … …  
    388203 END SUBROUTINE xios_read_field 
    389204 
    390  
    391   
    392  SUBROUTINE xios_write_field_scalar(name,field,nlev,nq) 
    393  USE genmod 
    394  USE mpipara 
    395  USE xios 
    396  USE grid_param 
    397  USE domain_mod 
    398  USE dimensions 
    399  USE spherical_geom_mod 
    400  USE geometry 
    401  USE mpi_mod 
    402  IMPLICIT NONE  
    403    CHARACTER(LEN=*),INTENT(IN) :: name 
    404    TYPE(t_field), POINTER :: field(:) 
    405    INTEGER,INTENT(IN) :: nlev 
    406    INTEGER,INTENT(IN) :: nq 
    407     
    408    REAL(rstd) :: field_tmp(ncell_i,nlev,nq) 
    409    TYPE(t_domain),POINTER :: d 
    410    INTEGER :: n,i,j,ij,ind 
    411     
    412    IF (field(1)%ndim==2) THEN 
    413      n=0 
    414      DO ind=1,ndomain 
    415         
    416        d=>domain(ind) 
    417          
    418        DO j=d%jj_begin,d%jj_end 
    419          DO i=d%ii_begin,d%ii_end 
    420            IF (d%own(i,j)) THEN 
    421              n=n+1 
    422              ij=d%iim*(j-1)+i 
    423              field_tmp(n,1,1)=field(ind)%rval2d(ij) 
    424            ENDIF 
    425          ENDDO 
    426        ENDDO 
    427      ENDDO 
    428    ELSE IF (field(1)%ndim==3) THEN 
    429      n=0 
    430      DO ind=1,ndomain 
    431        d=>domain(ind) 
    432          
    433        DO j=d%jj_begin,d%jj_end 
    434          DO i=d%ii_begin,d%ii_end 
    435            IF (d%own(i,j)) THEN 
    436              n=n+1 
    437              ij=d%iim*(j-1)+i 
    438              field_tmp(n,:,1)=field(ind)%rval3d(ij,:) 
    439            ENDIF 
    440          ENDDO 
    441        ENDDO 
    442      ENDDO 
    443    ELSE IF (field(1)%ndim==4) THEN 
    444      n=0 
    445      DO ind=1,ndomain 
    446        d=>domain(ind) 
    447          
    448        DO j=d%jj_begin,d%jj_end 
    449          DO i=d%ii_begin,d%ii_end 
    450            IF (d%own(i,j)) THEN 
    451              n=n+1 
    452              ij=d%iim*(j-1)+i 
    453              field_tmp(n,:,:)=field(ind)%rval4d(ij,:,:) 
    454            ENDIF 
    455          ENDDO 
    456        ENDDO 
    457      ENDDO      
    458    ENDIF 
    459     
    460    CALL xios_send_field(name,field_tmp) 
    461   
    462  END SUBROUTINE xios_write_field_scalar   
    463  
     205 SUBROUTINE xios_write_field_hex(name, field, cells, ncell_tot, nlev, nq) 
     206   CHARACTER(LEN=*),INTENT(IN) :: name 
     207   TYPE(t_field) :: field(:) 
     208   TYPE(t_cellset), TARGET :: cells(:) 
     209   INTEGER,INTENT(IN) :: ncell_tot, nlev, nq 
     210 
     211   REAL(rstd) :: field_tmp(ncell_tot,nlev,nq) 
     212   TYPE(t_cellset), POINTER :: cellset 
     213   INTEGER :: ind, n_beg, n_end, n, ij, sgn 
     214   LOGICAL :: signed 
     215 
     216   IF(ALLOCATED(cells(1)%sgn)) THEN 
     217      signed=.TRUE. 
     218   ELSE 
     219      signed=.FALSE. 
     220      sgn=1 
     221   END IF 
     222 
     223   n_beg=0 
     224   DO ind=1,ndomain 
     225      cellset => cells(ind) 
     226      n_end = n_beg + cellset%ncell 
     227      DO n=1, cellset%ncell 
     228         ij = cellset%ij(n) 
     229         IF(signed) sgn = cellset%sgn(n) 
     230         SELECT CASE(field(1)%ndim) 
     231         CASE(2) 
     232            field_tmp(n_beg+n,1,1) = sgn*field(ind)%rval2d(ij) 
     233         CASE(3) 
     234            field_tmp(n_beg+n,:,1) = sgn*field(ind)%rval3d(ij,:) 
     235         CASE(4) 
     236            field_tmp(n_beg+n,:,:) = sgn*field(ind)%rval4d(ij,:,:) 
     237         END SELECT 
     238      END DO 
     239   END DO 
     240   CALL xios_send_field(name,field_tmp)   
     241 END SUBROUTINE xios_write_field_hex 
     242 
     243SUBROUTINE xios_read_field_hex(name, field, cells, ncell_tot, nlev, nq) 
     244   CHARACTER(LEN=*),INTENT(IN) :: name 
     245   TYPE(t_field) :: field(:) 
     246   TYPE(t_cellset), TARGET :: cells(:) 
     247   INTEGER,INTENT(IN) :: ncell_tot, nlev, nq 
     248 
     249   REAL(rstd) :: field_tmp(ncell_tot,nlev,nq) 
     250   TYPE(t_cellset), POINTER :: cellset 
     251   INTEGER :: ind, n_beg, n_end, n, ij, sgn 
     252   LOGICAL :: signed 
     253 
     254   CALL xios_recv_field(name,field_tmp) 
     255 
     256   IF(ALLOCATED(cells(1)%sgn)) THEN 
     257      signed=.TRUE. 
     258   ELSE 
     259      signed=.FALSE. 
     260      sgn=1 
     261   END IF 
     262 
     263   n_beg=0 
     264   DO ind=1,ndomain 
     265      cellset => cells(ind) 
     266      n_end = n_beg + cellset%ncell 
     267      DO n=1, cellset%ncell 
     268         ij = cellset%ij(n) 
     269         IF(signed) sgn = cellset%sgn(n) 
     270         SELECT CASE(field(1)%ndim) 
     271         CASE(2) 
     272            field(ind)%rval2d(ij) = sgn*field_tmp(n_beg+n,1,1)  
     273         CASE(3) 
     274            field(ind)%rval3d(ij,:) = sgn*field_tmp(n_beg+n,:,1) 
     275         CASE(4) 
     276            field(ind)%rval4d(ij,:,:) = sgn*field_tmp(n_beg+n,:,:) 
     277         END SELECT 
     278      END DO 
     279   END DO 
     280 END SUBROUTINE xios_read_field_hex 
    464281 
    465282 SUBROUTINE xios_read_var(name,field) 
    … …  
    473290   CALL bcast_omp(field) 
    474291 END SUBROUTINE 
    475  
    476  SUBROUTINE xios_read_field_scalar(name,field,nlev,nq) 
    477  USE genmod 
    478  USE mpipara 
    479  USE xios 
    480  USE grid_param 
    481  USE domain_mod 
    482  USE dimensions 
    483  USE spherical_geom_mod 
    484  USE geometry 
    485  USE mpi_mod 
    486  IMPLICIT NONE  
    487    CHARACTER(LEN=*),INTENT(IN) :: name 
    488    TYPE(t_field), POINTER :: field(:) 
    489    INTEGER,INTENT(IN) :: nlev 
    490    INTEGER,INTENT(IN) :: nq 
    491     
    492    REAL(rstd) :: field_tmp(ncell_i,nlev,nq) 
    493    TYPE(t_domain),POINTER :: d 
    494    INTEGER :: n,i,j,ij,ind 
    495  
    496    CALL xios_recv_field(name,field_tmp) 
    497     
    498    IF (field(1)%ndim==2) THEN 
    499      n=0 
    500      DO ind=1,ndomain 
    501         
    502        d=>domain(ind) 
    503          
    504        DO j=d%jj_begin,d%jj_end 
    505          DO i=d%ii_begin,d%ii_end 
    506            IF (d%own(i,j)) THEN 
    507              n=n+1 
    508              ij=d%iim*(j-1)+i 
    509              field(ind)%rval2d(ij)=field_tmp(n,1,1) 
    510            ENDIF 
    511          ENDDO 
    512        ENDDO 
    513      ENDDO 
    514    ELSE IF (field(1)%ndim==3) THEN 
    515      n=0 
    516      DO ind=1,ndomain 
    517        d=>domain(ind) 
    518          
    519        DO j=d%jj_begin,d%jj_end 
    520          DO i=d%ii_begin,d%ii_end 
    521            IF (d%own(i,j)) THEN 
    522              n=n+1 
    523              ij=d%iim*(j-1)+i 
    524              field(ind)%rval3d(ij,:)=field_tmp(n,:,1) 
    525            ENDIF 
    526          ENDDO 
    527        ENDDO 
    528      ENDDO 
    529    ELSE IF (field(1)%ndim==4) THEN 
    530      n=0 
    531      DO ind=1,ndomain 
    532        d=>domain(ind) 
    533          
    534        DO j=d%jj_begin,d%jj_end 
    535          DO i=d%ii_begin,d%ii_end 
    536            IF (d%own(i,j)) THEN 
    537              n=n+1 
    538              ij=d%iim*(j-1)+i 
    539              field(ind)%rval4d(ij,:,:)=field_tmp(n,:,:) 
    540            ENDIF 
    541          ENDDO 
    542        ENDDO 
    543      ENDDO      
    544    ENDIF 
    545   
    546  END SUBROUTINE xios_read_field_scalar 
    547  
    548  SUBROUTINE xios_write_field_U(name,field,nlev,nq) 
    549  USE genmod 
    550  USE mpipara 
    551  USE xios 
    552  USE grid_param 
    553  USE domain_mod 
    554  USE dimensions 
    555  USE spherical_geom_mod 
    556  USE geometry 
    557  USE mpi_mod 
    558  IMPLICIT NONE  
    559    CHARACTER(LEN=*),INTENT(IN) :: name 
    560    TYPE(t_field), POINTER :: field(:) 
    561    INTEGER,INTENT(IN) :: nlev 
    562    INTEGER,INTENT(IN) :: nq 
    563     
    564    REAL(rstd) :: field_tmp(ncell_e,nlev,nq) 
    565    TYPE(t_domain),POINTER :: d 
    566    INTEGER :: n,i,j,k,ij,ind 
    567     
    568    IF (field(1)%ndim==2) THEN 
    569      n=0 
    570      DO ind=1,ndomain 
    571        d=>domain(ind) 
    572        CALL swap_dimensions(ind) 
    573        CALL swap_geometry(ind) 
    574  
    575        DO j=d%jj_begin,d%jj_end 
    576          DO i=d%ii_begin,d%ii_end 
    577            DO k=0,5 
    578              IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j & 
    579                                                                                  .AND. d%edge_assign_pos(k,i,j)==k) THEN 
    580                n=n+1 
    581                ij=iim*(j-1)+i 
    582                Field_tmp(n,1,1)=d%edge_assign_sign(k,i,j)*field(ind)%rval2d(ij+d%u_pos(k+1)) 
    583              ENDIF 
    584            ENDDO 
    585          ENDDO 
    586        ENDDO 
    587      ENDDO        
    588    
    589    ELSE IF (field(1)%ndim==3) THEN 
    590  
    591      n=0 
    592      DO ind=1,ndomain 
    593        d=>domain(ind) 
    594        CALL swap_dimensions(ind) 
    595        CALL swap_geometry(ind) 
    596  
    597        DO j=d%jj_begin,d%jj_end 
    598          DO i=d%ii_begin,d%ii_end 
    599            DO k=0,5 
    600              IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j & 
    601                                                                                   .AND. d%edge_assign_pos(k,i,j)==k) THEN 
    602                n=n+1 
    603                ij=iim*(j-1)+i 
    604                Field_tmp(n,:,1)=d%edge_assign_sign(k,i,j)*field(ind)%rval3d(ij+d%u_pos(k+1),:) 
    605              ENDIF 
    606            ENDDO 
    607          ENDDO 
    608        ENDDO 
    609      ENDDO        
    610  
    611    ELSE IF (field(1)%ndim==4) THEN 
    612  
    613      n=0 
    614      DO ind=1,ndomain 
    615        d=>domain(ind) 
    616        CALL swap_dimensions(ind) 
    617        CALL swap_geometry(ind) 
    618  
    619        DO j=d%jj_begin,d%jj_end 
    620          DO i=d%ii_begin,d%ii_end 
    621            DO k=0,5 
    622              IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j & 
    623                                                                                   .AND. d%edge_assign_pos(k,i,j)==k) THEN 
    624                n=n+1 
    625                ij=iim*(j-1)+i 
    626                Field_tmp(n,:,:)=d%edge_assign_sign(k,i,j)*field(ind)%rval4d(ij+d%u_pos(k+1),:,:) 
    627              ENDIF 
    628            ENDDO 
    629          ENDDO 
    630        ENDDO 
    631      ENDDO        
    632  
    633    ENDIF 
    634     
    635    CALL xios_send_field(name,field_tmp) 
    636   
    637  END SUBROUTINE xios_write_field_u  
    638  
    639  
    640  SUBROUTINE xios_read_field_u(name,field,nlev,nq) 
    641  USE genmod 
    642  USE mpipara 
    643  USE xios 
    644  USE grid_param 
    645  USE domain_mod 
    646  USE dimensions 
    647  USE spherical_geom_mod 
    648  USE geometry 
    649  USE mpi_mod 
    650  IMPLICIT NONE  
    651    CHARACTER(LEN=*),INTENT(IN) :: name 
    652    TYPE(t_field), POINTER :: field(:) 
    653    INTEGER,INTENT(IN) :: nlev 
    654    INTEGER,INTENT(IN) :: nq 
    655     
    656    REAL(rstd) :: field_tmp(ncell_e,nlev,nq) 
    657    TYPE(t_domain),POINTER :: d 
    658    INTEGER :: n,i,j,k,ij,ind 
    659  
    660    CALL xios_recv_field(name,field_tmp) 
    661     
    662    IF (field(1)%ndim==2) THEN 
    663      n=0 
    664      DO ind=1,ndomain 
    665        d=>domain(ind) 
    666        CALL swap_dimensions(ind) 
    667        CALL swap_geometry(ind) 
    668  
    669        DO j=d%jj_begin,d%jj_end 
    670          DO i=d%ii_begin,d%ii_end 
    671            DO k=0,5 
    672              IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j & 
    673                                                                                  .AND. d%edge_assign_pos(k,i,j)==k) THEN 
    674                n=n+1 
    675                ij=iim*(j-1)+i 
    676                field(ind)%rval2d(ij+d%u_pos(k+1))=Field_tmp(n,1,1)*d%edge_assign_sign(k,i,j) 
    677              ENDIF 
    678            ENDDO 
    679          ENDDO 
    680        ENDDO 
    681      ENDDO        
    682    
    683    ELSE IF (field(1)%ndim==3) THEN 
    684  
    685      n=0 
    686      DO ind=1,ndomain 
    687        d=>domain(ind) 
    688        CALL swap_dimensions(ind) 
    689        CALL swap_geometry(ind) 
    690  
    691        DO j=d%jj_begin,d%jj_end 
    692          DO i=d%ii_begin,d%ii_end 
    693            DO k=0,5 
    694              IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j & 
    695                                                                                   .AND. d%edge_assign_pos(k,i,j)==k) THEN 
    696                n=n+1 
    697                ij=iim*(j-1)+i 
    698                field(ind)%rval3d(ij+d%u_pos(k+1),:)=Field_tmp(n,:,1)*d%edge_assign_sign(k,i,j) 
    699              ENDIF 
    700            ENDDO 
    701          ENDDO 
    702        ENDDO 
    703      ENDDO        
    704  
    705    ELSE IF (field(1)%ndim==4) THEN 
    706  
    707      n=0 
    708      DO ind=1,ndomain 
    709        d=>domain(ind) 
    710        CALL swap_dimensions(ind) 
    711        CALL swap_geometry(ind) 
    712  
    713        DO j=d%jj_begin,d%jj_end 
    714          DO i=d%ii_begin,d%ii_end 
    715            DO k=0,5 
    716              IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j & 
    717                                                                                   .AND. d%edge_assign_pos(k,i,j)==k) THEN 
    718                n=n+1 
    719                ij=iim*(j-1)+i 
    720                field(ind)%rval4d(ij+d%u_pos(k+1),:,:)=Field_tmp(n,:,:)*d%edge_assign_sign(k,i,j) 
    721              ENDIF 
    722            ENDDO 
    723          ENDDO 
    724        ENDDO 
    725      ENDDO        
    726  
    727    ENDIF 
    728     
    729   
    730  END SUBROUTINE xios_read_field_u  
    731  
    732  
    733  
    734        
    735  SUBROUTINE xios_write_field_vort(name,field,nlev,iq) 
    736  USE genmod 
    737  USE mpipara 
    738  USE xios 
    739  USE grid_param 
    740  USE domain_mod 
    741  USE dimensions 
    742  USE spherical_geom_mod 
    743  USE geometry 
    744  USE mpi_mod 
    745  IMPLICIT NONE  
    746    CHARACTER(LEN=*),INTENT(IN) :: name 
    747    TYPE(t_field), POINTER :: field(:) 
    748    INTEGER,INTENT(IN) :: nlev 
    749    INTEGER,INTENT(IN),OPTIONAL :: iq 
    750     
    751    REAL(rstd) :: field_tmp(ncell_v,nlev) 
    752    TYPE(t_domain),POINTER :: d 
    753    INTEGER :: n,i,j,ij,ind 
    754     
    755    IF (field(1)%ndim==2) THEN 
    756      n=0 
    757      DO ind=1,ndomain 
    758        d=>domain(ind) 
    759        CALL swap_dimensions(ind)   
    760         
    761        DO j=d%jj_begin+1,d%jj_end 
    762          DO i=d%ii_begin,d%ii_end-1 
    763            n=n+1 
    764            ij=iim*(j-1)+i 
    765            Field_tmp(n,1)=field(ind)%rval2d(ij+z_down) 
    766          ENDDO 
    767        ENDDO 
    768  
    769        DO j=d%jj_begin,d%jj_end-1 
    770          DO i=d%ii_begin+1,d%ii_end 
    771            n=n+1 
    772            ij=iim*(j-1)+i 
    773            Field_tmp(n,1)=field(ind)%rval2d(ij+z_up) 
    774           ENDDO 
    775        ENDDO 
    776            
    777      ENDDO 
    778  
    779    ELSE IF (field(1)%ndim==3) THEN 
    780      n=0 
    781      DO ind=1,ndomain 
    782        d=>domain(ind) 
    783        CALL swap_dimensions(ind)     
    784               
    785        DO j=d%jj_begin+1,d%jj_end 
    786          DO i=d%ii_begin,d%ii_end-1 
    787            n=n+1 
    788            ij=iim*(j-1)+i 
    789            Field_tmp(n,:)=field(ind)%rval3d(ij+z_down,:) 
    790          ENDDO 
    791        ENDDO 
    792  
    793        DO j=d%jj_begin,d%jj_end-1 
    794          DO i=d%ii_begin+1,d%ii_end 
    795            n=n+1 
    796            ij=iim*(j-1)+i 
    797            Field_tmp(n,:)=field(ind)%rval3d(ij+z_up,:) 
    798           ENDDO 
    799        ENDDO 
    800            
    801      ENDDO 
    802  
    803    ELSE IF (field(1)%ndim==4) THEN 
    804      n=0 
    805      DO ind=1,ndomain 
    806        d=>domain(ind) 
    807        CALL swap_dimensions(ind)  
    808                 
    809        DO j=d%jj_begin+1,d%jj_end 
    810          DO i=d%ii_begin,d%ii_end-1 
    811            n=n+1 
    812            ij=iim*(j-1)+i 
    813            Field_tmp(n,:)=field(ind)%rval4d(ij+z_down,:,iq) 
    814          ENDDO 
    815        ENDDO 
    816  
    817        DO j=d%jj_begin,d%jj_end-1 
    818          DO i=d%ii_begin+1,d%ii_end 
    819            n=n+1 
    820            ij=iim*(j-1)+i 
    821            Field_tmp(n,:)=field(ind)%rval4d(ij+z_up,:,iq) 
    822           ENDDO 
    823        ENDDO 
    824            
    825      ENDDO 
    826  
    827    ENDIF 
    828     
    829    CALL xios_send_field(name,field_tmp) 
    830   
    831  END SUBROUTINE xios_write_field_vort  
    832  
    833  SUBROUTINE xios_read_field_vort(name,field,nlev,iq) 
    834  USE genmod 
    835  USE mpipara 
    836  USE xios 
    837  USE grid_param 
    838  USE domain_mod 
    839  USE dimensions 
    840  USE spherical_geom_mod 
    841  USE geometry 
    842  USE mpi_mod 
    843  IMPLICIT NONE  
    844    CHARACTER(LEN=*),INTENT(IN) :: name 
    845    TYPE(t_field), POINTER :: field(:) 
    846    INTEGER,INTENT(IN) :: nlev 
    847    INTEGER,INTENT(IN),OPTIONAL :: iq 
    848     
    849    REAL(rstd) :: field_tmp(ncell_v,nlev) 
    850    TYPE(t_domain),POINTER :: d 
    851    INTEGER :: n,i,j,ij,ind 
    852  
    853    CALL xios_recv_field(name,field_tmp) 
    854  
    855     
    856    IF (field(1)%ndim==2) THEN 
    857      n=0 
    858      DO ind=1,ndomain 
    859        d=>domain(ind) 
    860        CALL swap_dimensions(ind)   
    861         
    862        DO j=d%jj_begin+1,d%jj_end 
    863          DO i=d%ii_begin,d%ii_end-1 
    864            n=n+1 
    865            ij=iim*(j-1)+i 
    866            field(ind)%rval2d(ij+z_down)=Field_tmp(n,1) 
    867          ENDDO 
    868        ENDDO 
    869  
    870        DO j=d%jj_begin,d%jj_end-1 
    871          DO i=d%ii_begin+1,d%ii_end 
    872            n=n+1 
    873            ij=iim*(j-1)+i 
    874            Field_tmp(n,1)=field(ind)%rval2d(ij+z_up) 
    875            field(ind)%rval2d(ij+z_up)=Field_tmp(n,1) 
    876           ENDDO 
    877        ENDDO 
    878            
    879      ENDDO 
    880  
    881    ELSE IF (field(1)%ndim==3) THEN 
    882      n=0 
    883      DO ind=1,ndomain 
    884        d=>domain(ind) 
    885        CALL swap_dimensions(ind)     
    886               
    887        DO j=d%jj_begin+1,d%jj_end 
    888          DO i=d%ii_begin,d%ii_end-1 
    889            n=n+1 
    890            ij=iim*(j-1)+i 
    891            field(ind)%rval3d(ij+z_down,:)=Field_tmp(n,:) 
    892          ENDDO 
    893        ENDDO 
    894  
    895        DO j=d%jj_begin,d%jj_end-1 
    896          DO i=d%ii_begin+1,d%ii_end 
    897            n=n+1 
    898            ij=iim*(j-1)+i 
    899            field(ind)%rval3d(ij+z_up,:)=Field_tmp(n,:) 
    900           ENDDO 
    901        ENDDO 
    902            
    903      ENDDO 
    904  
    905    ELSE IF (field(1)%ndim==4) THEN 
    906      n=0 
    907      DO ind=1,ndomain 
    908        d=>domain(ind) 
    909        CALL swap_dimensions(ind)  
    910                 
    911        DO j=d%jj_begin+1,d%jj_end 
    912          DO i=d%ii_begin,d%ii_end-1 
    913            n=n+1 
    914            ij=iim*(j-1)+i 
    915            field(ind)%rval4d(ij+z_down,:,iq)=Field_tmp(n,:) 
    916          ENDDO 
    917        ENDDO 
    918  
    919        DO j=d%jj_begin,d%jj_end-1 
    920          DO i=d%ii_begin+1,d%ii_end 
    921            n=n+1 
    922            ij=iim*(j-1)+i 
    923            field(ind)%rval4d(ij+z_up,:,iq)=Field_tmp(n,:) 
    924           ENDDO 
    925        ENDDO 
    926            
    927      ENDDO 
    928  
    929    ENDIF 
    930   
    931  END SUBROUTINE xios_read_field_vort  
    932  
    933  
    934  
    935  
    936292   
    937293 SUBROUTINE xios_write_field_finalize 
    938  IMPLICIT NONE 
    939  
    940294!$OMP BARRIER 
    941295!$OMP MASTER  
    … …  
    943297!$OMP END MASTER 
    944298!$OMP BARRIER 
    945  
    946299 END SUBROUTINE xios_write_field_finalize 
    947300 
    948301 SUBROUTINE xios_set_context 
    949  IMPLICIT NONE    
    950302  TYPE(xios_context) :: ctx_hdl 
    951303 
    … …  
    971323   
    972324  SUBROUTINE xios_init 
    973    IMPLICIT NONE 
    974           
    975     using_xios=.FALSE. 
    976      
     325    using_xios=.FALSE.     
    977326  END SUBROUTINE xios_init 
    978327   
    979328  SUBROUTINE xios_send_field_scalar(name,field) 
    980   IMPLICIT NONE 
    981329    CHARACTER(LEN=*),INTENT(IN) :: name 
    982330    REAL,INTENT(IN) :: field 
    … …  
    984332 
    985333  SUBROUTINE xios_send_field_1d(name,field) 
    986   IMPLICIT NONE 
    987334    CHARACTER(LEN=*),INTENT(IN) :: name 
    988335    REAL,INTENT(IN) :: field(:) 
    … …  
    991338  SUBROUTINE xios_write_field(name,field) 
    992339  USE field_mod 
    993   IMPLICIT NONE  
    994340   CHARACTER(LEN=*),INTENT(IN) :: name 
    995341   TYPE(t_field), POINTER :: field(:) 
    … …  
    998344  SUBROUTINE xios_read_field(name,field) 
    999345  USE field_mod 
    1000   IMPLICIT NONE  
    1001346   CHARACTER(LEN=*),INTENT(IN) :: name 
    1002347   TYPE(t_field), POINTER :: field(:) 
    … …  
    1010355  
    1011356  SUBROUTINE xios_update_calendar(step) 
    1012   IMPLICIT NONE 
    1013357   INTEGER, INTENT(IN):: step   
    1014358  END SUBROUTINE xios_update_calendar 

  • codes/icosagcm/devel/src/sphere/compute_geometry.f90

    r863 r880  
    11MODULE compute_geometry_mod 
    22  USE geometry 
     3  USE dimensions 
     4 
     5  USE domain_mod, ONLY : t_domain, t_cellset, & 
     6       domain, ndomain, assigned_domain, & 
     7       domain_glo, ndomain_glo, domloc_glo_ind, swap_needed 
     8  USE omp_para, ONLY : is_omp_level_master, is_master 
     9  USE transfert_mod, ONLY : req_i0, req_i1, t_message, transfert_request, transfert_message, init_message 
     10 
     11  USE spherical_geom_mod, ONLY : xyz2lonlat, circumcenter, & 
     12       compute_centroid, centroid, & 
     13       surf_triangle, dist_cart, div_arc_bis, & 
     14       schmidt_transform 
     15  USE vector, ONLY : norm, cross_product2 
     16 
    317  IMPLICIT NONE 
    418  PRIVATE 
    … …  
    1024 
    1125  SUBROUTINE update_circumcenters  
    12     USE domain_mod 
    13     USE dimensions 
    14     USE spherical_geom_mod 
    15     USE vector 
    16     USE transfert_mod 
    17     USE omp_para 
    18  
    19     IMPLICIT NONE 
     26 
    2027    REAL(rstd) :: x1(3),x2(3) 
    2128    REAL(rstd) :: vect(3,6) 
    … …  
    3643     
    3744    DO ind=1,ndomain 
    38       IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
     45 
    3946      CALL swap_dimensions(ind) 
    4047      CALL swap_geometry(ind) 
    … …  
    5562 
    5663  SUBROUTINE remap_schmidt_loc 
    57     USE spherical_geom_mod 
    5864    USE getin_mod 
    59     USE omp_para 
    60     USE domain_mod 
    61     USE dimensions 
    62     IMPLICIT NONE 
    6365    INTEGER :: ind,i,j,n 
    6466    REAL(rstd) :: schmidt_factor, schmidt_lon, schmidt_lat 
    … …  
    9193 
    9294  SUBROUTINE optimize_geometry 
    93   USE metric 
    94   USE spherical_geom_mod 
    95   USE domain_mod 
    96   USE dimensions 
    97   USE transfert_mod 
    98   USE vector 
    9995  USE getin_mod 
    100   USE omp_para 
    101   IMPLICIT NONE 
    10296    INTEGER :: nb_it=0 
    10397    TYPE(t_domain),POINTER :: d 
    … …  
    193187    ! copy position of generators and vertices back into domain(:)%xyz/vertex 
    194188    ! so that XIOS/create_header_gen gets the right values 
    195     USE omp_para 
    196     USE dimensions 
    197     USE domain_mod  
    198     USE transfert_mpi_mod 
     189    USE transfert_mpi_mod, ONLY : gather_field, bcast_field 
    199190 
    200191    INTEGER :: ind,i,j,k,n 
    … …  
    253244  SUBROUTINE set_geometry 
    254245  USE metric 
    255   USE vector 
    256   USE spherical_geom_mod 
    257   USE domain_mod 
    258   USE dimensions 
    259   USE transfert_mod 
    260   USE getin_mod 
    261   USE omp_para 
    262   IMPLICIT NONE 
    263246 
    264247    REAL(rstd) :: surf(6) 
    … …  
    279262    CALL remap_schmidt_loc 
    280263    CALL update_circumcenters 
    281     ! copy position of generators and vertices back into domain(:)%xyz/vertex 
    282     ! so that XIOS gets the right values 
    283     CALL update_domain 
    284264 
    285265    DO ind=1,ndomain 
    … …  
    527507  
    528508  END SUBROUTINE set_geometry 
    529    
     509 
    530510  SUBROUTINE compute_wee(n,pos,w) 
    531   IMPLICIT NONE 
    532511    INTEGER,INTENT(IN) :: n 
    533512    INTEGER,INTENT(IN) :: pos 
    … …  
    555534   END SUBROUTINE compute_wee 
    556535             
    557  
    558536   
    559537  SUBROUTINE compute_geometry 
    560538    USE grid_param 
    561     USE domain_mod, ONLY : swap_needed 
    562539    USE init_unstructured_mod, ONLY : read_local_mesh 
    563     IMPLICIT NONE 
    564  
     540    USE set_bounds_mod, ONLY : set_bounds 
    565541    CALL allocate_geometry 
    566542 
    … …  
    568544    CASE(grid_ico) 
    569545       CALL set_geometry 
     546       ! copy position of generators and vertices back into domain_glo(:)%xyz/vertex 
     547       ! so that write_field gets the right values 
     548       CALL update_domain 
     549       CALL set_bounds(domain_glo, .TRUE.) 
     550       CALL set_bounds(domain, .FALSE.) 
    570551    CASE(grid_unst) 
    571552       swap_needed=.FALSE. 
Note: See TracChangeset for help on using the changeset viewer.