Changeset 26


Ignore:
Timestamp:
07/26/12 15:25:40 (13 years ago)
Author:
ymipsl
Message:

Implementation of parallelism
implementation of iterative laplacian for dissipation

Location:
codes/icosagcm/trunk/src
Files:
5 added
1 deleted
8 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/caldyn_gcm.f90

    r22 r26  
    401401    
    402402!!!  Compute mass flux 
    403 !! question à thomas : meilleure pondération de la masse sur les liens ? 
     403!! question ï¿œ thomas : meilleure pondï¿œration de la masse sur les liens ? 
    404404 
    405405  DO l = 1, llm 

  • codes/icosagcm/trunk/src/dissip_gcm.f90

    r22 r26  
    22  USE icosa 
    33 
    4   TYPE(t_field),POINTER,SAVE :: f_gradrot(:) 
    5   TYPE(t_request),POINTER :: req_dissip(:)  
    6   TYPE(t_field),POINTER,SAVE :: f_due(:) 
     4  PRIVATE 
     5 
     6  TYPE(t_field),POINTER,SAVE :: f_due_diss1(:) 
     7  TYPE(t_field),POINTER,SAVE :: f_due_diss2(:) 
     8 
     9  TYPE(t_field),POINTER,SAVE :: f_theta(:) 
     10  TYPE(t_field),POINTER,SAVE :: f_dtheta_diss(:) 
     11  TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz_diss(:) 
     12   
    713  
    8   INTEGER,PARAMETER :: nitergdiv=1 
    9   INTEGER,PARAMETER :: nitergrot=1 
    10   INTEGER,PARAMETER :: niterdivgrad=1 
     14  INTEGER,SAVE :: nitergdiv=1 
     15  INTEGER,SAVE :: nitergrot=1 
     16  INTEGER,SAVE :: niterdivgrad=1 
    1117 
    1218  REAL,ALLOCATABLE,SAVE :: tau_graddiv(:) 
    … …  
    1824  REAL,SAVE :: cdivgrad 
    1925   
    20   INTEGER :: idissip 
    21   REAL    :: dtdissip 
    22    
     26  INTEGER,SAVE :: idissip 
     27  REAL,SAVE    :: dtdissip 
     28   
     29  PUBLIC init_dissip, dissip 
    2330CONTAINS 
    2431 
    … …  
    2633  USE icosa 
    2734  IMPLICIT NONE   
    28     CALL allocate_field(f_gradrot,field_u,type_real,llm) 
    29     CALL allocate_field(f_due,field_u,type_real,llm) 
     35    CALL allocate_field(f_due_diss1,field_u,type_real,llm) 
     36    CALL allocate_field(f_due_diss2,field_u,type_real,llm) 
     37    CALL allocate_field(f_theta,field_t,type_real,llm) 
     38    CALL allocate_field(f_dtheta_diss,field_t,type_real,llm) 
     39    CALL allocate_field(f_dtheta_rhodz_diss,field_t,type_real,llm) 
     40     
    3041    ALLOCATE(tau_graddiv(llm)) 
    3142    ALLOCATE(tau_gradrot(llm))     
    … …  
    3647  USE icosa 
    3748  USE disvert_mod 
     49  USE mpi_mod 
     50  USE mpipara 
    3851   
    3952  IMPLICIT NONE 
    … …  
    5669   
    5770             
    58   INTEGER :: i,j,n,ind,it 
     71  INTEGER :: i,j,n,ind,it,iter 
    5972 
    6073    CALL allocate_dissip 
    … …  
    6780    CALL getin("tau_graddiv",tau) 
    6881    tau_graddiv(:)=tau 
     82 
     83    CALL getin("nitergdiv",nitergdiv) 
    6984     
    7085    tau_gradrot(:)=5000 
    71     CALL getin("tau_gradrot",tau) 
     86    CALL getin("tau_gradrot",tau_gradrot) 
    7287    tau_gradrot(:)=tau 
     88 
     89    CALL getin("nitergrot",nitergrot) 
    7390     
    7491 
    … …  
    7693    CALL getin("tau_divgrad",tau) 
    7794    tau_divgrad(:)=tau 
    78    
    79     CALL create_request(field_u,req_dissip) 
    80  
    81     DO ind=1,ndomain 
    82       DO i=ii_begin,ii_end 
    83         CALL request_add_point(ind,i,jj_begin-1,req_dissip,rup) 
    84         CALL request_add_point(ind,i+1,jj_begin-1,req_dissip,lup) 
    85       ENDDO 
    86  
    87       DO j=jj_begin,jj_end 
    88         CALL request_add_point(ind,ii_end+1,j,req_dissip,left) 
    89         CALL request_add_point(ind,ii_end+1,j-1,req_dissip,lup) 
    90       ENDDO     
    91      
    92       DO i=ii_begin,ii_end 
    93         CALL request_add_point(ind,i,jj_end+1,req_dissip,ldown) 
    94         CALL request_add_point(ind,i-1,jj_end+1,req_dissip,rdown) 
    95       ENDDO     
    96  
    97       DO j=jj_begin,jj_end 
    98         CALL request_add_point(ind,ii_begin-1,j,req_dissip,right) 
    99         CALL request_add_point(ind,ii_begin-1,j+1,req_dissip,rdown) 
    100       ENDDO    
    101  
    102       DO i=ii_begin+1,ii_end-1 
    103         CALL request_add_point(ind,i,jj_begin,req_dissip,right) 
    104         CALL request_add_point(ind,i,jj_end,req_dissip,right) 
    105       ENDDO 
    106      
    107       DO j=jj_begin+1,jj_end-1 
    108         CALL request_add_point(ind,ii_begin,j,req_dissip,rup) 
    109         CALL request_add_point(ind,ii_end,j,req_dissip,rup) 
    110       ENDDO    
    111  
    112       CALL request_add_point(ind,ii_begin+1,jj_begin,req_dissip,left) 
    113       CALL request_add_point(ind,ii_begin,jj_begin+1,req_dissip,ldown) 
    114       CALL request_add_point(ind,ii_begin+1,jj_end,req_dissip,left) 
    115       CALL request_add_point(ind,ii_end,jj_begin+1,req_dissip,ldown) 
    116       
    117     ENDDO 
     95 
     96    CALL getin("niterdivgrad",niterdivgrad) 
     97   
     98!   CALL create_request(field_u,req_dissip) 
     99 
     100!   DO ind=1,ndomain 
     101!     DO i=ii_begin,ii_end 
     102!       CALL request_add_point(ind,i,jj_begin-1,req_dissip,rup) 
     103!       CALL request_add_point(ind,i+1,jj_begin-1,req_dissip,lup) 
     104!     ENDDO 
     105 
     106!     DO j=jj_begin,jj_end 
     107!       CALL request_add_point(ind,ii_end+1,j,req_dissip,left) 
     108!       CALL request_add_point(ind,ii_end+1,j-1,req_dissip,lup) 
     109!     ENDDO     
     110!    
     111!     DO i=ii_begin,ii_end 
     112!       CALL request_add_point(ind,i,jj_end+1,req_dissip,ldown) 
     113!       CALL request_add_point(ind,i-1,jj_end+1,req_dissip,rdown) 
     114!     ENDDO     
     115 
     116!     DO j=jj_begin,jj_end 
     117!       CALL request_add_point(ind,ii_begin-1,j,req_dissip,right) 
     118!       CALL request_add_point(ind,ii_begin-1,j+1,req_dissip,rdown) 
     119!     ENDDO    
     120 
     121!     DO i=ii_begin+1,ii_end-1 
     122!       CALL request_add_point(ind,i,jj_begin,req_dissip,right) 
     123!       CALL request_add_point(ind,i,jj_end,req_dissip,right) 
     124!     ENDDO 
     125!    
     126!     DO j=jj_begin+1,jj_end-1 
     127!       CALL request_add_point(ind,ii_begin,j,req_dissip,rup) 
     128!       CALL request_add_point(ind,ii_end,j,req_dissip,rup) 
     129!     ENDDO    
     130 
     131!     CALL request_add_point(ind,ii_begin+1,jj_begin,req_dissip,left) 
     132!     CALL request_add_point(ind,ii_begin,jj_begin+1,req_dissip,ldown) 
     133!     CALL request_add_point(ind,ii_begin+1,jj_end,req_dissip,left) 
     134!     CALL request_add_point(ind,ii_end,jj_begin+1,req_dissip,ldown) 
     135!     
     136!   ENDDO 
    118137 
    119138 
    … …  
    145164 
    146165 
    147     DO it=1,20 
    148  
    149       CALL transfert_request(f_u,req_dissip) 
    150       CALL transfert_request(f_u,req_dissip) 
    151  
     166    DO it=1,500 
     167       
    152168      dumax=0 
    153       DO ind=1,ndomain 
    154         CALL swap_dimensions(ind) 
    155         CALL swap_geometry(ind) 
    156         u=f_u(ind) 
    157         du=f_du(ind) 
    158         CALL gradiv(u,du,1) 
    159       ENDDO 
    160       CALL transfert_request(f_du,req_dissip) 
    161       CALL transfert_request(f_du,req_dissip) 
     169      DO iter=1,nitergdiv 
     170        CALL transfert_request(f_u,req_e1) 
     171        DO ind=1,ndomain 
     172          CALL swap_dimensions(ind) 
     173          CALL swap_geometry(ind) 
     174          u=f_u(ind) 
     175          du=f_du(ind) 
     176          CALL compute_gradiv(u,du,1) 
     177          u=du 
     178        ENDDO 
     179      ENDDO 
     180 
     181      CALL transfert_request(f_du,req_e1) 
    162182       
    163183      DO ind=1,ndomain 
    … …  
    166186        u=f_u(ind) 
    167187        du=f_du(ind) 
    168          
     188          
    169189        DO j=jj_begin,jj_end 
    170190          DO i=ii_begin,ii_end 
    … …  
    177197      ENDDO 
    178198 
     199      IF (using_mpi) CALL MPI_ALLREDUCE(dumax,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 
     200 
    179201      DO ind=1,ndomain 
    180202        CALL swap_dimensions(ind) 
    … …  
    212234 
    213235 
    214     DO it=1,20 
     236    DO it=1,500 
    215237  
    216       CALL transfert_request(f_u,req_dissip) 
    217       CALL transfert_request(f_u,req_dissip) 
    218  
    219238      dumax=0 
    220       DO ind=1,ndomain 
    221         CALL swap_dimensions(ind) 
    222         CALL swap_geometry(ind) 
    223         u=f_u(ind) 
    224         du=f_du(ind) 
    225         CALL gradrot(u,du,1) 
    226       ENDDO 
    227  
    228       CALL transfert_request(f_du,req_dissip) 
    229       CALL transfert_request(f_du,req_dissip) 
     239      DO iter=1,nitergrot 
     240        CALL transfert_request(f_u,req_e1) 
     241        DO ind=1,ndomain 
     242          CALL swap_dimensions(ind) 
     243          CALL swap_geometry(ind) 
     244          u=f_u(ind) 
     245          du=f_du(ind) 
     246          CALL compute_gradrot(u,du,1) 
     247          u=du 
     248        ENDDO 
     249      ENDDO 
     250 
     251      CALL transfert_request(f_du,req_e1) 
    230252       
    231253      DO ind=1,ndomain 
    … …  
    244266        ENDDO 
    245267      ENDDO 
     268 
     269      IF (using_mpi) CALL MPI_ALLREDUCE(dumax,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 
    246270       
    247271      DO ind=1,ndomain 
    … …  
    253277      ENDDO 
    254278       
    255 !      PRINT *,"gradrot : it :",it ,": dumax",(dumax+dumaxm1)/2 
    256279      PRINT *,"gradrot : it :",it ,": dumax",dumax 
    257280 
    … …  
    259282    PRINT *,"gradrot : dumax",dumax 
    260283   
    261     cgradrot=dumax**(-1/nitergrot)  
     284    cgradrot=dumax**(-1./nitergrot)  
    262285    PRINT *,"cgradrot : ",cgradrot 
    263286    
    … …  
    279302    ENDDO 
    280303 
    281     DO it=1,20 
    282  
    283       CALL transfert_request(f_theta,req_i1) 
     304    DO it=1,500 
    284305 
    285306      dthetamax=0 
    286       DO ind=1,ndomain 
    287         CALL swap_dimensions(ind) 
    288         CALL swap_geometry(ind) 
    289         theta=f_theta(ind) 
    290         dtheta=f_dtheta(ind) 
    291         CALL divgrad(theta,dtheta,1) 
     307      DO iter=1,niterdivgrad 
     308        CALL transfert_request(f_theta,req_i1) 
     309        DO ind=1,ndomain 
     310          CALL swap_dimensions(ind) 
     311          CALL swap_geometry(ind) 
     312          theta=f_theta(ind) 
     313          dtheta=f_dtheta(ind) 
     314          CALL compute_divgrad(theta,dtheta,1) 
     315          theta=dtheta 
     316        ENDDO 
    292317      ENDDO 
    293318!      CALL writefield("divgrad",f_dtheta) 
    … …  
    308333        ENDDO 
    309334      ENDDO 
    310  
     335      IF (using_mpi) CALL MPI_ALLREDUCE(dthetamax,dthetamax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 
    311336      PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax 
    312337 
    … …  
    323348    PRINT *,"divgrad : divgrad",dthetamax 
    324349   
    325     cdivgrad=dthetamax**(-1/niterdivgrad)  
     350    cdivgrad=dthetamax**(-1./niterdivgrad)  
    326351    PRINT *,"cdivgrad : ",cdivgrad 
    327352     
    … …  
    358383    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    359384    TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 
    360     REAL(rstd),POINTER         :: ue(:,:) 
     385 
    361386    REAL(rstd),POINTER         :: due(:,:) 
    362     REAL(rstd),POINTER         :: ps(:) 
    363     REAL(rstd),POINTER         :: theta_rhodz(:,:) 
     387    REAL(rstd),POINTER         :: due_diss1(:,:) 
     388    REAL(rstd),POINTER         :: due_diss2(:,:) 
    364389    REAL(rstd),POINTER         :: dtheta_rhodz(:,:) 
    365  
    366     REAL(rstd)                 :: theta(iim*jjm,llm) 
    367     REAL(rstd)                 :: du_dissip(3*iim*jjm,llm) 
    368     REAL(rstd)                 :: dtheta_dissip(iim*jjm,llm) 
    369     REAL(rstd)                 :: dtheta_rhodz_dissip(iim*jjm,llm) 
     390    REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:) 
    370391 
    371392    INTEGER :: ind 
    372393    INTEGER :: l,i,j,n 
    373394     
    374     CALL transfert_request(f_ue,req_e1) 
    375     CALL transfert_request(f_ue,req_e1) 
    376     CALL transfert_request(f_theta_rhodz,req_i1) 
    377     CALL transfert_request(f_ps,req_i1) 
    378  
     395    CALL gradiv(f_ue,f_due_diss1) 
     396    CALL gradrot(f_ue,f_due_diss2) 
     397    CALL theta_rhodz2theta(f_ps,f_theta_rhodz,f_theta) 
     398    CALL divgrad(f_theta,f_dtheta_diss) 
     399    CALL theta2theta_rhodz(f_ps,f_dtheta_diss,f_dtheta_rhodz_diss) 
     400     
     401    DO ind=1,ndomain 
     402      CALL swap_dimensions(ind) 
     403      CALL swap_geometry(ind) 
     404      due=f_due(ind)  
     405      due_diss1=f_due_diss1(ind) 
     406      due_diss2=f_due_diss2(ind) 
     407      dtheta_rhodz=f_dtheta_rhodz(ind) 
     408      dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind) 
     409             
     410      DO l=1,llm 
     411        DO j=jj_begin,jj_end 
     412          DO i=ii_begin,ii_end 
     413            n=(j-1)*iim+i 
     414 
     415            due(n+u_right,l)=due(n+u_right,l)-0.5*(due_diss1(n+u_right,l)/tau_graddiv(l)+due_diss2(n+u_right,l)/tau_gradrot(l)) 
     416            due(n+u_lup,l)=due(n+u_lup,l)    -0.5*(due_diss1(n+u_lup,l)  /tau_graddiv(l)+due_diss2(n+u_lup,l)  /tau_gradrot(l)) 
     417            due(n+u_ldown,l)=due(n+u_ldown,l)-0.5*(due_diss1(n+u_ldown,l)/tau_graddiv(l)+due_diss2(n+u_ldown,l)/tau_gradrot(l)) 
     418 
     419            dtheta_rhodz(n,l)=dtheta_rhodz(n,l)-0.5*dtheta_rhodz_diss(n,l)/tau_divgrad(l) 
     420 
     421          ENDDO 
     422        ENDDO 
     423      ENDDO 
     424    ENDDO 
     425       
     426  END SUBROUTINE dissip 
     427 
     428  SUBROUTINE gradiv(f_ue,f_due) 
     429  USE icosa 
     430  IMPLICIT NONE 
     431    TYPE(t_field),POINTER :: f_ue(:) 
     432    TYPE(t_field),POINTER :: f_due(:) 
     433    REAL(rstd),POINTER    :: ue(:,:) 
     434    REAL(rstd),POINTER    :: due(:,:) 
     435    INTEGER :: ind 
     436    INTEGER :: it 
     437        
    379438    DO ind=1,ndomain 
    380439      CALL swap_dimensions(ind) 
    … …  
    382441      ue=f_ue(ind) 
    383442      due=f_due(ind)  
    384       ps=f_ps(ind)  
    385       theta_rhodz=f_theta_rhodz(ind) 
    386       dtheta_rhodz=f_dtheta_rhodz(ind)  
    387  
    388 !$OMP PARALLEL DEFAULT(SHARED) 
    389      CALL compute_dissip(ue,due,ps,theta_rhodz,dtheta_rhodz) 
    390 !$OMP END PARALLEL 
    391        
    392     ENDDO 
    393        
    394   END SUBROUTINE dissip 
    395  
    396  
    397   SUBROUTINE compute_dissip(ue,due,ps,theta_rhodz,dtheta_rhodz) 
    398   USE icosa 
    399   USE theta2theta_rhodz_mod 
     443      due=ue 
     444    ENDDO 
     445 
     446    DO it=1,nitergdiv 
     447         
     448      CALL transfert_request(f_due,req_e1) 
     449         
     450      DO ind=1,ndomain 
     451        CALL swap_dimensions(ind) 
     452        CALL swap_geometry(ind) 
     453        due=f_due(ind)  
     454        CALL compute_gradiv(due,due,llm) 
     455      ENDDO 
     456    ENDDO 
     457 
     458  END SUBROUTINE gradiv 
     459   
     460 
     461  SUBROUTINE gradrot(f_ue,f_due) 
     462  USE icosa 
    400463  IMPLICIT NONE 
    401     REAL(rstd)         :: ue(3*iim*jjm,llm) 
    402     REAL(rstd)         :: due(3*iim*jjm,llm) 
    403     REAL(rstd)         :: ps(iim*jjm) 
    404     REAL(rstd)         :: theta_rhodz(iim*jjm,llm) 
    405     REAL(rstd)         :: dtheta_rhodz(iim*jjm,llm) 
    406  
    407     REAL(rstd),SAVE,ALLOCATABLE :: theta(:,:) 
    408     REAL(rstd),SAVE,ALLOCATABLE :: du_dissip(:,:) 
    409     REAL(rstd),SAVE,ALLOCATABLE :: dtheta_dissip(:,:) 
    410     REAL(rstd),SAVE,ALLOCATABLE :: dtheta_rhodz_dissip(:,:) 
    411  
     464    TYPE(t_field),POINTER :: f_ue(:) 
     465    TYPE(t_field),POINTER :: f_due(:) 
     466    REAL(rstd),POINTER    :: ue(:,:) 
     467    REAL(rstd),POINTER    :: due(:,:) 
    412468    INTEGER :: ind 
    413     INTEGER :: l,i,j,n 
    414  
    415 !$OMP BARRIER 
    416 !$OMP MASTER 
    417       ALLOCATE(theta(iim*jjm,llm)) 
    418       ALLOCATE(du_dissip(3*iim*jjm,llm)) 
    419       ALLOCATE(dtheta_dissip(iim*jjm,llm)) 
    420       ALLOCATE(dtheta_rhodz_dissip(iim*jjm,llm)) 
    421 !$OMP END MASTER 
    422 !$OMP BARRIER 
    423      
    424       CALL gradiv(ue,du_dissip,llm) 
    425       DO l=1,llm 
    426 !$OMP DO 
    427         DO j=jj_begin,jj_end 
    428           DO i=ii_begin,ii_end 
    429             n=(j-1)*iim+i 
    430             due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_graddiv(l)*0.5 
    431             due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_graddiv(l)*0.5 
    432             due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_graddiv(l)*0.5 
    433           ENDDO 
    434         ENDDO 
    435       ENDDO 
    436        
    437       CALL gradrot(ue,du_dissip,llm) 
    438  
    439       DO l=1,llm 
    440 !$OMP DO 
    441         DO j=jj_begin,jj_end 
    442           DO i=ii_begin,ii_end 
    443             n=(j-1)*iim+i 
    444             due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_gradrot(l)*0.5 
    445             due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_gradrot(l)*0.5 
    446             due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_gradrot(l)*0.5 
    447           ENDDO 
    448         ENDDO 
    449       ENDDO  
    450        
    451       CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1) 
    452       CALL divgrad(theta,dtheta_dissip,llm) 
    453       CALL compute_theta2theta_rhodz(ps,dtheta_dissip,dtheta_rhodz_dissip,0) 
    454  
    455       DO l=1,llm 
    456 !$OMP DO 
    457         DO j=jj_begin,jj_end 
    458           DO i=ii_begin,ii_end 
    459             n=(j-1)*iim+i 
    460             dtheta_rhodz(n,l)=dtheta_rhodz(n,l)+dtheta_rhodz_dissip(n,l)/tau_divgrad(l)*0.5 
    461           ENDDO 
    462         ENDDO 
    463       ENDDO  
    464  
    465 !$OMP BARRIER 
    466 !$OMP MASTER 
    467       DEALLOCATE(theta) 
    468       DEALLOCATE(du_dissip) 
    469       DEALLOCATE(dtheta_dissip) 
    470       DEALLOCATE(dtheta_rhodz_dissip) 
    471 !$OMP END MASTER 
    472 !$OMP BARRIER       
    473  
    474   END SUBROUTINE compute_dissip 
    475    
    476    
    477   SUBROUTINE gradiv(ue,gradivu_e,ll) 
     469    INTEGER :: it 
     470        
     471    DO ind=1,ndomain 
     472      CALL swap_dimensions(ind) 
     473      CALL swap_geometry(ind) 
     474      ue=f_ue(ind) 
     475      due=f_due(ind)  
     476      due=ue 
     477    ENDDO 
     478 
     479    DO it=1,nitergrot 
     480         
     481      CALL transfert_request(f_due,req_e1) 
     482         
     483      DO ind=1,ndomain 
     484        CALL swap_dimensions(ind) 
     485        CALL swap_geometry(ind) 
     486        due=f_due(ind)  
     487        CALL compute_gradrot(due,due,llm) 
     488      ENDDO 
     489 
     490    ENDDO 
     491 
     492  END SUBROUTINE gradrot 
     493   
     494  SUBROUTINE divgrad(f_theta,f_dtheta) 
     495  USE icosa 
     496  IMPLICIT NONE 
     497    TYPE(t_field),POINTER :: f_theta(:) 
     498    TYPE(t_field),POINTER :: f_dtheta(:) 
     499    REAL(rstd),POINTER    :: theta(:,:) 
     500    REAL(rstd),POINTER    :: dtheta(:,:) 
     501    INTEGER :: ind 
     502    INTEGER :: it 
     503        
     504    DO ind=1,ndomain 
     505      CALL swap_dimensions(ind) 
     506      CALL swap_geometry(ind) 
     507      theta=f_theta(ind) 
     508      dtheta=f_dtheta(ind)  
     509      dtheta=theta 
     510    ENDDO 
     511 
     512    DO it=1,niterdivgrad 
     513         
     514      CALL transfert_request(f_dtheta,req_i1) 
     515         
     516      DO ind=1,ndomain 
     517        CALL swap_dimensions(ind) 
     518        CALL swap_geometry(ind) 
     519        dtheta=f_dtheta(ind)  
     520        CALL compute_divgrad(dtheta,dtheta,llm) 
     521      ENDDO 
     522 
     523    ENDDO 
     524 
     525  END SUBROUTINE divgrad 
     526    
     527   
     528 
     529!  SUBROUTINE compute_dissip(ue,due,ps,theta_rhodz,dtheta_rhodz) 
     530!  USE icosa 
     531!  USE theta2theta_rhodz_mod 
     532!  IMPLICIT NONE 
     533!    REAL(rstd)         :: ue(3*iim*jjm,llm) 
     534!    REAL(rstd)         :: due(3*iim*jjm,llm) 
     535!    REAL(rstd)         :: ps(iim*jjm) 
     536!    REAL(rstd)         :: theta_rhodz(iim*jjm,llm) 
     537!    REAL(rstd)         :: dtheta_rhodz(iim*jjm,llm) 
     538! 
     539!    REAL(rstd),SAVE,ALLOCATABLE :: theta(:,:) 
     540!    REAL(rstd),SAVE,ALLOCATABLE :: du_dissip(:,:) 
     541!    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_dissip(:,:) 
     542!    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_rhodz_dissip(:,:) 
     543! 
     544!    INTEGER :: ind 
     545!    INTEGER :: l,i,j,n 
     546! 
     547!!$OMP BARRIER 
     548!!$OMP MASTER 
     549!      ALLOCATE(theta(iim*jjm,llm)) 
     550!      ALLOCATE(du_dissip(3*iim*jjm,llm)) 
     551!      ALLOCATE(dtheta_dissip(iim*jjm,llm)) 
     552!      ALLOCATE(dtheta_rhodz_dissip(iim*jjm,llm)) 
     553!!$OMP END MASTER 
     554!!$OMP BARRIER 
     555!     
     556!      CALL gradiv(ue,du_dissip,llm) 
     557!      DO l=1,llm 
     558!!$OMP DO 
     559!        DO j=jj_begin,jj_end 
     560!          DO i=ii_begin,ii_end 
     561!            n=(j-1)*iim+i 
     562!            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_graddiv(l)*0.5 
     563!            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_graddiv(l)*0.5 
     564!            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_graddiv(l)*0.5 
     565!          ENDDO 
     566!        ENDDO 
     567!      ENDDO 
     568!       
     569!      CALL gradrot(ue,du_dissip,llm) 
     570! 
     571!      DO l=1,llm 
     572!!$OMP DO 
     573!        DO j=jj_begin,jj_end 
     574!          DO i=ii_begin,ii_end 
     575!            n=(j-1)*iim+i 
     576!            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_gradrot(l)*0.5 
     577!            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_gradrot(l)*0.5 
     578!            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_gradrot(l)*0.5 
     579!          ENDDO 
     580!        ENDDO 
     581!      ENDDO  
     582!       
     583!      CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1) 
     584!      CALL divgrad(theta,dtheta_dissip,llm) 
     585!      CALL compute_theta2theta_rhodz(ps,dtheta_dissip,dtheta_rhodz_dissip,0) 
     586! 
     587!      DO l=1,llm 
     588!!$OMP DO 
     589!        DO j=jj_begin,jj_end 
     590!          DO i=ii_begin,ii_end 
     591!            n=(j-1)*iim+i 
     592!            dtheta_rhodz(n,l)=dtheta_rhodz(n,l)+dtheta_rhodz_dissip(n,l)/tau_divgrad(l)*0.5 
     593!          ENDDO 
     594!        ENDDO 
     595!      ENDDO  
     596! 
     597!!$OMP BARRIER 
     598!!$OMP MASTER 
     599!      DEALLOCATE(theta) 
     600!      DEALLOCATE(du_dissip) 
     601!      DEALLOCATE(dtheta_dissip) 
     602!      DEALLOCATE(dtheta_rhodz_dissip) 
     603!!$OMP END MASTER 
     604!!$OMP BARRIER       
     605! 
     606!  END SUBROUTINE compute_dissip 
     607   
     608   
     609  SUBROUTINE compute_gradiv(ue,gradivu_e,ll) 
    478610  USE icosa 
    479611  IMPLICIT NONE 
    … …  
    528660        DO i=ii_begin,ii_end 
    529661          n=(j-1)*iim+i 
    530           gradivu_e(n+u_right,l)=gradivu_e(n+u_right,l)*cgraddiv 
    531           gradivu_e(n+u_lup,l)=gradivu_e(n+u_lup,l)*cgraddiv 
    532           gradivu_e(n+u_ldown,l)=gradivu_e(n+u_ldown,l)*cgraddiv 
     662          gradivu_e(n+u_right,l)=-gradivu_e(n+u_right,l)*cgraddiv 
     663          gradivu_e(n+u_lup,l)=-gradivu_e(n+u_lup,l)*cgraddiv 
     664          gradivu_e(n+u_ldown,l)=-gradivu_e(n+u_ldown,l)*cgraddiv 
    533665        ENDDO 
    534666      ENDDO 
    … …  
    542674 
    543675     
    544   END SUBROUTINE gradiv 
    545    
    546   SUBROUTINE divgrad(theta,divgrad_i,ll) 
     676  END SUBROUTINE compute_gradiv 
     677   
     678  SUBROUTINE compute_divgrad(theta,divgrad_i,ll) 
    547679  USE icosa 
    548680  IMPLICIT NONE 
    … …  
    598730        DO i=ii_begin,ii_end 
    599731          n=(j-1)*iim+i 
    600           divgrad_i(n,l)=divgrad_i(n,l)*cdivgrad 
     732          divgrad_i(n,l)=-divgrad_i(n,l)*cdivgrad 
    601733        ENDDO 
    602734      ENDDO 
    … …  
    609741!$OMP BARRIER 
    610742 
    611   END SUBROUTINE divgrad 
    612  
    613      
    614   SUBROUTINE gradrot(ue,gradrot_e,ll) 
     743  END SUBROUTINE compute_divgrad 
     744 
     745     
     746  SUBROUTINE compute_gradrot(ue,gradrot_e,ll) 
    615747  USE icosa 
    616748  IMPLICIT NONE 
    … …  
    668800          n=(j-1)*iim+i 
    669801     
    670           gradrot_e(n+u_right,l)=gradrot_e(n+u_right,l)*cgradrot        
    671           gradrot_e(n+u_lup,l)=gradrot_e(n+u_lup,l)*cgradrot 
    672           gradrot_e(n+u_ldown,l)=gradrot_e(n+u_ldown,l)*cgradrot 
     802          gradrot_e(n+u_right,l)=-gradrot_e(n+u_right,l)*cgradrot        
     803          gradrot_e(n+u_lup,l)=-gradrot_e(n+u_lup,l)*cgradrot 
     804          gradrot_e(n+u_ldown,l)=-gradrot_e(n+u_ldown,l)*cgradrot 
    673805         
    674806        ENDDO 
    … …  
    682814!$OMP BARRIER  
    683815    
    684   END SUBROUTINE gradrot 
     816  END SUBROUTINE compute_gradrot 
    685817 
    686818 

  • codes/icosagcm/trunk/src/domain.f90

    r21 r26  
    5959  TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain(:) 
    6060  INTEGER :: ndomain 
    61  
    62  
     61  TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain_glo(:) 
     62  INTEGER :: ndomain_glo 
     63 
     64  INTEGER,ALLOCATABLE,SAVE :: domglo_rank(:) 
     65  INTEGER,ALLOCATABLE,SAVE :: domglo_loc_ind(:) 
     66  INTEGER,ALLOCATABLE,SAVE :: domloc_glo_ind(:) 
     67   
    6368CONTAINS 
    6469 
    … …  
    7075  TYPE(t_domain),POINTER :: d 
    7176  
    72     ndomain=nsplit_i*nsplit_j*nb_face 
    73     ALLOCATE(domain(ndomain)) 
    74    
     77    ndomain_glo=nsplit_i*nsplit_j*nb_face 
     78    ALLOCATE(domain_glo(ndomain_glo)) 
     79    ALLOCATE(domglo_rank(ndomain_glo)) 
     80    ALLOCATE(domglo_loc_ind(ndomain_glo)) 
     81 
    7582    ind=0 
    7683    DO nf=1,nb_face 
    … …  
    7885        DO ni=1,nsplit_i 
    7986          ind=ind+1 
    80           d=>domain(ind) 
     87          d=>domain_glo(ind) 
    8188          quotient=(iim_glo/nsplit_i) 
    8289          rest=MOD(iim_glo,nsplit_i) 
    … …  
    139146  END SUBROUTINE create_domain 
    140147   
     148  SUBROUTINE copy_domain(d1,d2) 
     149  IMPLICIT NONE 
     150  INTEGER :: face 
     151  TYPE(t_domain),TARGET,INTENT(IN) :: d1 
     152  TYPE(t_domain), INTENT(OUT) :: d2 
     153   
     154    d2%iim = d1%iim 
     155    d2%jjm = d1%jjm 
     156    d2%ii_begin = d1%ii_begin 
     157    d2%jj_begin = d1%jj_begin 
     158    d2%ii_end = d1%ii_end 
     159    d2%jj_end = d1%jj_end 
     160    d2%ii_nb =  d1%ii_nb 
     161    d2%jj_nb = d1%jj_nb 
     162    d2%ii_begin_glo = d1%ii_begin_glo 
     163    d2%jj_begin_glo = d1%jj_begin_glo 
     164    d2%ii_end_glo = d1%ii_end_glo 
     165    d2%jj_end_glo = d1%jj_end_glo 
     166    d2%assign_domain => d1%assign_domain 
     167    d2%assign_i => d1%assign_i 
     168    d2%assign_j => d1%assign_j 
     169    d2%edge_assign_domain => d1%edge_assign_domain 
     170    d2%edge_assign_i => d1%edge_assign_i 
     171    d2%edge_assign_j => d1%edge_assign_j 
     172    d2%edge_assign_pos => d1%edge_assign_pos 
     173    d2%xyz => d1%xyz 
     174    d2%neighbour => d1%neighbour 
     175    d2%delta => d1%delta 
     176    d2%vertex => d1%vertex 
     177    d2%own => d1%own 
     178    d2%ne => d1%ne 
     179     
     180    d2%t_right = d1%t_right 
     181    d2%t_rup = d1%t_rup 
     182    d2%t_lup = d1%t_lup 
     183    d2%t_left = d1%t_left 
     184    d2%t_ldown = d1%t_ldown 
     185    d2%t_rdown = d1%t_rdown 
     186 
     187    d2%u_right = d1%u_right 
     188    d2%u_rup = d1%u_rup 
     189    d2%u_lup = d1%u_lup 
     190    d2%u_left = d1%u_left 
     191    d2%u_ldown = d1%u_ldown 
     192    d2%u_rdown = d1%u_rdown 
     193 
     194    d2%z_rup = d1%z_rup 
     195    d2%z_up = d1%z_up 
     196    d2%z_lup = d1%z_lup 
     197    d2%z_ldown = d1%z_ldown 
     198    d2%z_down = d1%z_down 
     199    d2%z_rdown = d1%z_rdown 
     200    
     201    d2%t_pos = d1%t_pos 
     202    d2%u_pos = d1%u_pos 
     203    d2%z_pos = d1%z_pos 
     204     
     205  END SUBROUTINE copy_domain 
     206   
    141207   
    142208  SUBROUTINE assign_cell 
    … …  
    150216      
    151217     
    152     DO ind_d=1,ndomain 
    153       d=>domain(ind_d) 
     218    DO ind_d=1,ndomain_glo 
     219      d=>domain_glo(ind_d) 
    154220      nf=d%face 
    155221      DO j=d%jj_begin,d%jj_end 
    … …  
    175241     
    176242     
    177     DO ind_d=1,ndomain 
    178       d=>domain(ind_d) 
     243    DO ind_d=1,ndomain_glo 
     244      d=>domain_glo(ind_d) 
    179245      nf=d%face 
    180246      DO j=d%jj_begin-1,d%jj_end+1 
    … …  
    231297    REAL(rstd) :: ng1(3),ng2(3)  
    232298 
    233     DO ind_d=1,ndomain 
    234       d=>domain(ind_d) 
     299    DO ind_d=1,ndomain_glo 
     300      d=>domain_glo(ind_d) 
    235301      DO j=d%jj_begin-1,d%jj_end+1 
    236302        DO i=d%ii_begin-1,d%ii_end+1 
    … …  
    252318    TYPE(t_domain),POINTER :: d  
    253319     
    254     DO ind_d=1,ndomain 
    255       d=>domain(ind_d) 
     320    DO ind_d=1,ndomain_glo 
     321      d=>domain_glo(ind_d) 
    256322      d%t_right=1 
    257323      d%t_left=-1 
    … …  
    301367  END SUBROUTINE set_neighbour_indice 
    302368       
    303        
     369  SUBROUTINE assign_domain 
     370  USE mpipara 
     371  IMPLICIT NONE 
     372    INTEGER :: nb_domain(0:mpi_size-1) 
     373    INTEGER :: rank, ind,ind_glo 
     374     
     375    DO rank=0,mpi_size-1 
     376      nb_domain(rank)=ndomain_glo/mpi_size 
     377      IF ( rank < MOD(ndomain_glo,mpi_size) ) nb_domain(rank)=nb_domain(rank)+1 
     378    ENDDO 
     379     
     380    ndomain=nb_domain(mpi_rank) 
     381    ALLOCATE(domain(ndomain)) 
     382    ALLOCATE(domloc_glo_ind(ndomain)) 
     383     
     384    rank=0 
     385    ind=0 
     386    DO ind_glo=1,ndomain_glo 
     387      ind=ind+1 
     388      domglo_rank(ind_glo)=rank 
     389      domglo_loc_ind(ind_glo)=ind 
     390      IF (rank==mpi_rank) THEN  
     391        CALL copy_domain(domain_glo(ind_glo),domain(ind)) 
     392        domloc_glo_ind(ind)=ind_glo 
     393      ENDIF 
     394       
     395      IF (ind==nb_domain(rank)) THEN 
     396        rank=rank+1 
     397        ind=0 
     398      ENDIF 
     399    ENDDO 
     400     
     401  END SUBROUTINE  assign_domain       
     402     
    304403           
    305404  SUBROUTINE compute_domain 
    306405  IMPLICIT NONE 
    307    
     406    CALL init_domain_param 
    308407    CALL create_domain 
    309408    CALL assign_cell 
    310409    CALL compute_boundary 
    311410    CALL set_neighbour_indice 
     411    CALL assign_domain 
    312412       
    313413  END SUBROUTINE compute_domain 

  • codes/icosagcm/trunk/src/domain_param.f90

    r12 r26  
    11MODULE domain_param 
    22 
    3 INTEGER, PARAMETER :: nsplit_i=1 
    4 INTEGER, PARAMETER :: nsplit_j=1 
    5 INTEGER, PARAMETER :: halo=2 
     3INTEGER :: nsplit_i 
     4INTEGER :: nsplit_j 
     5INTEGER :: halo=2 
    66 
     7INTEGER, PARAMETER :: default_nsplit_i=1 
     8INTEGER, PARAMETER :: default_nsplit_j=1 
    79 
     10CONTAINS 
     11   
     12  SUBROUTINE init_domain_param 
     13  USE ioipsl 
     14  IMPLICIT NONE 
     15    nsplit_i=default_nsplit_i 
     16    nsplit_j=default_nsplit_j 
     17    CALL getin('nsplit_i',nsplit_i) 
     18    CALL getin('nsplit_j',nsplit_j) 
     19  END SUBROUTINE init_domain_param 
     20   
    821END MODULE domain_param 

  • codes/icosagcm/trunk/src/field.f90

    r12 r26  
    2626    INTEGER :: field_type 
    2727    INTEGER :: data_type  
     28    INTEGER :: dim3 
     29    INTEGER :: dim4 
    2830  END TYPE t_field    
    2931 
    … …  
    5961      IF (PRESENT(dim2)) THEN 
    6062        field(ind)%ndim=4  
     63        field(ind)%dim4=dim2  
    6164      ELSE IF (PRESENT(dim1)) THEN 
    6265        field(ind)%ndim=3 
     66        field(ind)%dim3=dim1 
    6367      ELSE 
    6468        field(ind)%ndim=2 
    … …  
    96100    
    97101  END SUBROUTINE allocate_field 
    98    
     102 
     103  SUBROUTINE allocate_field_glo(field,field_type,data_type,dim1,dim2) 
     104  USE domain_mod 
     105  IMPLICIT NONE 
     106    TYPE(t_field),POINTER :: field(:) 
     107    INTEGER,INTENT(IN) :: field_type 
     108    INTEGER,INTENT(IN) :: data_type 
     109    INTEGER,OPTIONAL :: dim1,dim2 
     110    INTEGER :: ind 
     111    INTEGER :: ii_size,jj_size 
     112 
     113    ALLOCATE(field(ndomain_glo))     
     114 
     115    DO ind=1,ndomain_glo 
     116   
     117      IF (PRESENT(dim2)) THEN 
     118        field(ind)%ndim=4  
     119        field(ind)%dim4=dim2  
     120      ELSE IF (PRESENT(dim1)) THEN 
     121        field(ind)%ndim=3 
     122        field(ind)%dim3=dim1  
     123      ELSE 
     124        field(ind)%ndim=2 
     125      ENDIF 
     126     
     127     
     128      field(ind)%data_type=data_type 
     129      field(ind)%field_type=field_type 
     130     
     131      IF (field_type==field_T) THEN  
     132        jj_size=domain_glo(ind)%jjm 
     133      ELSE IF (field_type==field_U) THEN  
     134        jj_size=3*domain_glo(ind)%jjm 
     135      ELSE IF (field_type==field_Z) THEN  
     136        jj_size=2*domain_glo(ind)%jjm 
     137      ENDIF 
     138       
     139      ii_size=domain_glo(ind)%iim 
     140         
     141      IF (field(ind)%ndim==4) THEN 
     142        IF (data_type==type_integer) ALLOCATE(field(ind)%ival4d(ii_size*jj_size,dim1,dim2)) 
     143        IF (data_type==type_real)    ALLOCATE(field(ind)%rval4d(ii_size*jj_size,dim1,dim2)) 
     144        IF (data_type==type_logical) ALLOCATE(field(ind)%lval4d(ii_size*jj_size,dim1,dim2)) 
     145      ELSE IF (field(ind)%ndim==3) THEN 
     146        IF (data_type==type_integer) ALLOCATE(field(ind)%ival3d(ii_size*jj_size,dim1)) 
     147        IF (data_type==type_real)    ALLOCATE(field(ind)%rval3d(ii_size*jj_size,dim1)) 
     148        IF (data_type==type_logical) ALLOCATE(field(ind)%lval3d(ii_size*jj_size,dim1)) 
     149      ELSE IF (field(ind)%ndim==2) THEN 
     150        IF (data_type==type_integer) ALLOCATE(field(ind)%ival2d(ii_size*jj_size)) 
     151        IF (data_type==type_real)    ALLOCATE(field(ind)%rval2d(ii_size*jj_size)) 
     152        IF (data_type==type_logical) ALLOCATE(field(ind)%lval2d(ii_size*jj_size)) 
     153      ENDIF 
     154       
     155   ENDDO 
     156   
     157  END SUBROUTINE allocate_field_glo 
     158 
     159  SUBROUTINE deallocate_field(field) 
     160  USE domain_mod 
     161  IMPLICIT NONE 
     162    TYPE(t_field),POINTER :: field(:) 
     163    INTEGER :: data_type 
     164    INTEGER :: ind 
     165 
     166    DO ind=1,ndomain 
     167 
     168      data_type=field(ind)%data_type 
     169         
     170      IF (field(ind)%ndim==4) THEN 
     171        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival4d) 
     172        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval4d) 
     173        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval4d) 
     174      ELSE IF (field(ind)%ndim==3) THEN 
     175        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival3d) 
     176        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval3d) 
     177        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval3d) 
     178      ELSE IF (field(ind)%ndim==2) THEN 
     179        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival2d) 
     180        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval2d) 
     181        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval2d) 
     182      ENDIF 
     183       
     184   ENDDO 
     185   DEALLOCATE(field) 
     186        
     187  END SUBROUTINE deallocate_field 
     188 
     189  SUBROUTINE deallocate_field_glo(field) 
     190  USE domain_mod 
     191  IMPLICIT NONE 
     192    TYPE(t_field),POINTER :: field(:) 
     193    INTEGER :: data_type 
     194    INTEGER :: ind 
     195 
     196    DO ind=1,ndomain_glo 
     197 
     198      data_type=field(ind)%data_type 
     199         
     200      IF (field(ind)%ndim==4) THEN 
     201        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival4d) 
     202        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval4d) 
     203        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval4d) 
     204      ELSE IF (field(ind)%ndim==3) THEN 
     205        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival3d) 
     206        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval3d) 
     207        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval3d) 
     208      ELSE IF (field(ind)%ndim==2) THEN 
     209        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival2d) 
     210        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval2d) 
     211        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval2d) 
     212      ENDIF 
     213       
     214   ENDDO 
     215   DEALLOCATE(field) 
     216        
     217  END SUBROUTINE deallocate_field_glo 
     218     
    99219  SUBROUTINE getval_r2d(field_pt,field) 
    100220  IMPLICIT NONE   

  • codes/icosagcm/trunk/src/icosa_gcm.f90

    r19 r26  
    22  USE icosa 
    33  USE timeloop_gcm_mod 
    4   USE transfert_mod 
    54  USE disvert_mod 
    65  USE etat0_mod 
    76  USE wind_mod 
     7  USE mpipara 
    88  IMPLICIT NONE 
    99   
    1010  TYPE(t_field),POINTER :: sum_ne(:) 
     11  TYPE(t_field),POINTER :: sum_ne_glo(:) 
    1112  REAL(rstd),POINTER :: pt_sum_ne(:) 
    1213   
    … …  
    1617  REAL(rstd) :: centr(3),dist 
    1718  
     19  CALL init_mpipara 
     20   
    1821  CALL init_grid_param 
    1922  CALL compute_metric 
    2023  CALL compute_domain 
    2124  CALL init_transfert 
     25!  CALL allocate_field(sum_ne,field_T,type_real) 
     26!  CALL allocate_field_glo(sum_ne_glo,field_T,type_real) 
     27!   
     28! DO ind=1,ndomain 
     29!   CALL swap_dimensions(ind) 
     30!   pt_sum_ne=sum_ne(ind) 
     31!   DO j=jj_begin,jj_end 
     32!     DO i=ii_begin,ii_end     
     33!       n=(j-1)*iim+i 
     34!       pt_sum_ne(n)=domloc_glo_ind(ind) 
     35!     ENDDO 
     36!   ENDDO 
     37! ENDDO 
     38 
     39! CALL WriteField("domain",sum_ne) 
     40! CALL WriteField_mpi("domain",sum_ne) 
     41! CALL transfert_request(sum_ne,req_i1) 
     42! CALL WriteField_mpi("domain",sum_ne) 
     43! CALL close_files 
     44! CALL finalize_mpipara 
     45! STOP 
     46   
    2247  CALL compute_geometry 
    2348  CALL init_disvert  
    … …  
    5277  CALL WriteField("Ai",geom%Ai) 
    5378!  CALL WriteField("sum_ne",sum_ne) 
    54    
    5579  CALL timeloop 
     80 
     81  CALL close_files 
     82  CALL finalize_mpipara 
    5683   
    5784END PROGRAM ICOSA_gcm  

  • codes/icosagcm/trunk/src/icosa_sw.f90

    r19 r26  
    22  USE icosa 
    33  USE timeloop_sw_mod 
    4   USE transfert_mod 
    54  USE dissip_sw_mod 
    65  USE disvert_mod 

  • codes/icosagcm/trunk/src/write_field.f90

    r21 r26  
    66    INTEGER :: size 
    77    INTEGER,POINTER :: nc_id(:) 
     8    INTEGER :: displ 
    89  END TYPE ncvar 
    910 
    … …  
    3738      enddo 
    3839    end function GetFieldIndex 
    39   
    40  
    41        
    42     subroutine WriteField_gen(name,Field,dimx,dimy,dimz) 
    43     implicit none 
    44 !    include 'netcdf.inc' 
    45       character(len=*) :: name 
    46       integer :: dimx,dimy,dimz 
    47       real,dimension(dimx,dimy,dimz) :: Field 
    48       integer,dimension(dimx*dimy*dimz) :: ndex 
    49       integer :: status 
    50       integer :: index 
    51       integer :: start(4) 
    52       integer :: count(4) 
    53        
    54             
     40     
     41    SUBROUTINE Writefield(name_in,field,nind) 
     42    USE domain_mod 
     43    USE field_mod 
     44    USE transfert_mpi_mod 
     45    USE dimensions 
     46    USE mpipara 
     47    IMPLICIT NONE   
     48     CHARACTER(LEN=*),INTENT(IN) :: name_in 
     49      TYPE(t_field),POINTER :: field(:) 
     50      INTEGER,OPTIONAL,INTENT(IN) :: nind 
     51      TYPE(t_field),POINTER :: field_glo(:) 
     52       
     53      IF (field(1)%ndim==2) THEN 
     54        CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type) 
     55      ELSE IF (field(1)%ndim==3) THEN 
     56        CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3)  
     57      ELSE IF (field(1)%ndim==4) THEN 
     58        CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3,field(1)%dim4) 
     59      ENDIF 
     60       
     61      CALL gather_field(field,field_glo) 
     62       
     63      IF (mpi_rank==0) THEN 
     64        IF (PRESENT(nind)) THEN 
     65         CALL writefield_gen(name_in,field_glo,domain_glo,nind) 
     66        ELSE 
     67         CALL writefield_gen(name_in,field_glo,domain_glo,1,ndomain_glo) 
     68        ENDIF 
     69      ENDIF 
     70       
     71      CALL deallocate_field(field_glo) 
     72       
     73   END SUBROUTINE writefield 
     74    
     75!   SUBROUTINE Writefield(name_in,field,nind) 
     76!   USE netcdf 
     77!   USE domain_mod 
     78!   use field_mod 
     79!   USE dimensions 
     80!   USE geometry 
     81!   IMPLICIT NONE 
     82!     CHARACTER(LEN=*),INTENT(IN) :: name_in 
     83!     TYPE(t_field),POINTER :: field(:) 
     84!     INTEGER,OPTIONAL,INTENT(IN) :: nind 
     85!     REAL(r8),ALLOCATABLE :: field_val2d(:) 
     86!     REAL(r8),ALLOCATABLE :: field_val3d(:,:) 
     87!     REAL(r8),ALLOCATABLE :: field_val4d(:,:,:) 
     88!     TYPE(t_domain),POINTER :: d 
     89!     INTEGER :: Index 
     90!     INTEGER :: ind,i,j,k,n,ncell,q 
     91!     INTEGER :: iie,jje,iin,jjn 
     92!     INTEGER :: status 
     93!     CHARACTER(len=255) :: name 
     94!     CHARACTER(len=255) :: str_ind 
     95!     INTEGER :: ind_b,ind_e 
     96!     INTEGER :: halo_size 
     97!     LOGICAL :: single 
     98!      
     99!      
     100!     name=TRIM(ADJUSTL(name_in)) 
     101 
     102!     IF (PRESENT(nind)) THEN 
     103!       name=TRIM(name)//"_"//TRIM(int2str(nind)) 
     104!       PRINT *,"NAME",nind,int2str(nind),name 
     105!       ind_b=nind 
     106!       ind_e=nind 
     107!       halo_size=1 
     108!       single=.TRUE. 
     109!     ELSE 
     110!       ind_b=1 
     111!       ind_e=ndomain 
     112!       halo_size=0 
     113!       single=.FALSE. 
     114!     ENDIF       
     115 
     116!     Index=GetFieldIndex(name) 
     117!     if (Index==-1) then 
     118!       call create_header(name,field,nind) 
     119!       Index=GetFieldIndex(name) 
     120!     else 
     121!       FieldIndex(Index)=FieldIndex(Index)+1. 
     122!     endif 
     123!      
     124!     IF (Field(ind_b)%field_type==field_T) THEN 
     125!       ncell=1 
     126!       DO ind=ind_b,ind_e 
     127!         d=>domain(ind) 
     128!         IF (Field(ind)%field_type/=field_T) THEN 
     129!           PRINT *,"Writefield, grille non geree" 
     130!           RETURN 
     131!         ENDIF 
     132 
     133!       n=0 
     134!       DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     135!         DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     136!           IF (d%own(i,j) .OR. single) n=n+1 
     137!         ENDDO 
     138!       ENDDO 
     139 
     140!       IF (field(ind)%ndim==2) THEN 
     141!         ALLOCATE(Field_val2d(n)) 
     142!         n=0 
     143!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     144!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     145!             k=d%iim*(j-1)+i 
     146!             IF (d%own(i,j) .OR. single) THEN 
     147!               n=n+1 
     148!               Field_val2d(n)=field(ind)%rval2d(k) 
     149!             ENDIF 
     150!           ENDDO 
     151!         ENDDO 
     152!         status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  & 
     153!                             start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
     154!         DEALLOCATE(field_val2d) 
     155!       ELSE IF (field(ind)%ndim==3) THEN 
     156!         ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2))) 
     157!         n=0 
     158!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     159!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     160!             k=d%iim*(j-1)+i 
     161!             IF (d%own(i,j) .OR. single) THEN 
     162!               n=n+1 
     163!               Field_val3d(n,:)=field(ind)%rval3d(k,:) 
     164!             ENDIF 
     165!           ENDDO 
     166!         ENDDO 
     167!          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
     168!                              count=(/n,size(field(1)%rval3d,2),1 /)) 
     169!         DEALLOCATE(field_val3d) 
     170!       ELSE IF (field(1)%ndim==4) THEN 
     171 
     172!         DO q=1,FieldVarId(index)%size 
     173!          
     174!           ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2))) 
     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!               k=d%iim*(j-1)+i 
     179!               IF (d%own(i,j) .OR. single) THEN 
     180!                 n=n+1 
     181!                 Field_val3d(n,:)=field(ind)%rval4d(k,:,q) 
     182!               ENDIF 
     183!             ENDDO 
     184!           ENDDO 
     185!           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
     186!                              count=(/n,size(field(1)%rval4d,2),1 /)) 
     187!           DEALLOCATE(field_val3d) 
     188!         ENDDO          
     189!       ENDIF 
     190!        
     191!        PRINT *,NF90_STRERROR(status) 
     192!       ncell=ncell+n 
     193 
     194!    ENDDO 
     195!     
     196!    ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
     197!       ncell=1 
     198!       n=0 
     199!       DO ind=ind_b,ind_e 
     200!         d=>domain(ind) 
     201!         CALL swap_geometry(ind) 
     202!         CALL swap_dimensions(ind) 
     203!  
     204!         n=0 
     205!         DO j=jj_begin+1,jj_end 
     206!           DO i=ii_begin,ii_end-1 
     207!             n=n+1 
     208!           ENDDO 
     209!         ENDDO 
     210 
     211!         DO j=jj_begin,jj_end-1 
     212!           DO i=ii_begin+1,ii_end 
     213!             n=n+1 
     214!           ENDDO 
     215!         ENDDO 
     216 
     217!       IF (field(ind)%ndim==2) THEN 
     218!         ALLOCATE(Field_val2d(n)) 
     219 
     220!         n=0 
     221!         DO j=jj_begin+1,jj_end 
     222!           DO i=ii_begin,ii_end-1 
     223!             n=n+1 
     224!             k=iim*(j-1)+i 
     225!             Field_val2d(n)=field(ind)%rval2d(k+z_down) 
     226!           ENDDO 
     227!         ENDDO 
     228 
     229!         DO j=jj_begin,jj_end-1 
     230!           DO i=ii_begin+1,ii_end 
     231!             n=n+1 
     232!             k=iim*(j-1)+i 
     233!             Field_val2d(n)=field(ind)%rval2d(k+z_up) 
     234!           ENDDO 
     235!         ENDDO 
     236 
     237!         status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       & 
     238!                             Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
     239!         DEALLOCATE(field_val2d) 
     240 
     241!       ELSE IF (field(ind)%ndim==3) THEN 
     242!         ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2))) 
     243!         n=0 
     244!         DO j=jj_begin+1,jj_end 
     245!           DO i=ii_begin,ii_end-1 
     246!             n=n+1 
     247!             k=iim*(j-1)+i 
     248!             Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:) 
     249!           ENDDO 
     250!         ENDDO 
     251 
     252!         DO j=jj_begin,jj_end-1 
     253!           DO i=ii_begin+1,ii_end 
     254!             n=n+1 
     255!             k=iim*(j-1)+i 
     256!             Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:) 
     257!           ENDDO 
     258!         ENDDO 
     259!          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
     260!                              count=(/n,size(field(1)%rval3d,2),1 /)) 
     261!         DEALLOCATE(field_val3d) 
     262!       ELSE IF (field(1)%ndim==4) THEN 
     263 
     264!         DO q=1,FieldVarId(index)%size 
     265!           ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2))) 
     266!           n=0 
     267!           DO j=jj_begin+1,jj_end 
     268!             DO i=ii_begin,ii_end-1 
     269!               n=n+1 
     270!               k=iim*(j-1)+i 
     271!               Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q) 
     272!             ENDDO 
     273!           ENDDO 
     274 
     275!           DO j=jj_begin,jj_end-1 
     276!             DO i=ii_begin+1,ii_end 
     277!               n=n+1 
     278!               k=iim*(j-1)+i 
     279!               Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q) 
     280!             ENDDO 
     281!           ENDDO 
     282 
     283!           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,1,FieldIndex(Index) /),   & 
     284!                               count=(/n,size(field(1)%rval4d,2),1 /)) 
     285!           DEALLOCATE(field_val3d) 
     286!         ENDDO 
     287!       ENDIF 
     288!        
     289!        PRINT *,NF90_STRERROR(status) 
     290!       ncell=ncell+n 
     291 
     292!    ENDDO 
     293!     
     294!    ENDIF 
     295!    status=NF90_SYNC(FieldId(Index)) 
     296!      
     297!   END SUBROUTINE Writefield 
     298 
     299 
     300    SUBROUTINE Writefield_gen(name_in, field, domain_type, ind_b_in, ind_e_in ) 
     301    USE netcdf_mod 
     302    USE domain_mod 
     303    USE field_mod 
     304!    USE dimensions 
     305!    USE geometry 
     306    IMPLICIT NONE 
     307      CHARACTER(LEN=*),INTENT(IN) :: name_in 
     308      TYPE(t_field), POINTER :: field(:) 
     309      TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:) 
     310      INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in 
     311      INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in 
     312      REAL(r8),ALLOCATABLE :: field_val2d(:) 
     313      REAL(r8),ALLOCATABLE :: field_val3d(:,:) 
     314      REAL(r8),ALLOCATABLE :: field_val4d(:,:,:) 
     315      TYPE(t_domain),POINTER :: d 
     316      INTEGER :: Index 
     317      INTEGER :: ind,i,j,k,n,ncell,q 
     318      INTEGER :: iie,jje,iin,jjn 
     319      INTEGER :: status 
     320      CHARACTER(len=255) :: name 
     321      CHARACTER(len=255) :: str_ind 
     322      INTEGER :: ind_b,ind_e 
     323      INTEGER :: halo_size 
     324      LOGICAL :: single 
     325       
     326       
     327      name=TRIM(ADJUSTL(name_in)) 
     328 
     329      IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN 
     330        name=TRIM(name)//"_"//TRIM(int2str(ind_b)) 
     331        ind_b=ind_b_in 
     332        ind_e=ind_b_in 
     333        halo_size=1 
     334        single=.TRUE. 
     335      ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
     336        name=TRIM(name)//"_"//TRIM(int2str(ind_e)) 
     337        ind_b=ind_e_in 
     338        ind_e=ind_e_in 
     339        halo_size=1 
     340        single=.TRUE. 
     341      ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
     342        ind_b=ind_b_in 
     343        ind_e=ind_e_in 
     344        halo_size=0 
     345        single=.FALSE. 
     346      ELSE 
     347        halo_size=0 
     348        ind_b=1 
     349        ind_e=ndomain 
     350        single=.FALSE. 
     351      ENDIF       
     352       
    55353      Index=GetFieldIndex(name) 
    56354      if (Index==-1) then 
    57 !        call CreateNewField(name,dimx,dimy,dimz) 
    58         Index=GetFieldIndex(name) 
     355        call create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in) 
     356        Index=GetFieldIndex(name) 
    59357      else 
    60358        FieldIndex(Index)=FieldIndex(Index)+1. 
    61359      endif 
    62360       
    63       start(1)=1 
    64       start(2)=1 
    65       start(3)=1 
    66       start(4)=FieldIndex(Index) 
    67  
    68       count(1)=dimx 
    69       count(2)=dimy 
    70       count(3)=dimz 
    71       count(4)=1 
    72  
    73 !      status = NF_PUT_VARA_DOUBLE(FieldId(Index),FieldVarId(Index),start,count,Field) 
    74 !      status = NF_SYNC(FieldId(Index)) 
    75        
    76     end subroutine WriteField_gen 
    77  
    78  
    79     SUBROUTINE Writefield(name_in,field,nind) 
    80     USE netcdf 
     361      IF (Field(ind_b)%field_type==field_T) THEN 
     362 
     363        ncell=0 
     364        DO ind=ind_b,ind_e 
     365          d=>domain_type(ind) 
     366          DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     367            DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     368              IF (d%assign_domain(i,j)==ind .OR. single) ncell=ncell+1 
     369            ENDDO 
     370          ENDDO 
     371        ENDDO 
     372         
     373        IF (field(1)%ndim==2) THEN 
     374          ALLOCATE(Field_val2d(ncell)) 
     375          n=0 
     376          DO ind=ind_b,ind_e 
     377            d=>domain_type(ind) 
     378            DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     379              DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     380                k=d%iim*(j-1)+i 
     381                IF (d%assign_domain(i,j)==ind .OR. single) THEN 
     382                  n=n+1 
     383                  Field_val2d(n)=field(ind)%rval2d(k) 
     384                ENDIF 
     385              ENDDO 
     386            ENDDO 
     387          ENDDO 
     388          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  & 
     389                              start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /)) 
     390          DEALLOCATE(field_val2d) 
     391        ELSE IF (field(1)%ndim==3) THEN 
     392          ALLOCATE(Field_val3d(ncell,size(field(1)%rval3d,2))) 
     393          n=0 
     394          DO ind=ind_b,ind_e 
     395            d=>domain_type(ind) 
     396            DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     397              DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     398                k=d%iim*(j-1)+i 
     399                IF (d%assign_domain(i,j)==ind .OR. single) THEN 
     400                  n=n+1 
     401                  Field_val3d(n,:)=field(ind)%rval3d(k,:) 
     402                ENDIF 
     403              ENDDO 
     404            ENDDO 
     405          ENDDO 
     406            
     407          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1,FieldIndex(Index) /),   & 
     408                               count=(/ncell,size(field(1)%rval3d,2),1 /)) 
     409          DEALLOCATE(field_val3d) 
     410        ELSE IF (field(1)%ndim==4) THEN 
     411 
     412          DO q=1,FieldVarId(index)%size 
     413           
     414            ALLOCATE(Field_val3d(n,size(field(1)%rval4d,2))) 
     415            n=0 
     416            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_val3d(n,:)=field(ind)%rval4d(k,:,q) 
     424                  ENDIF 
     425                ENDDO 
     426              ENDDO 
     427            ENDDO 
     428             
     429            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,FieldIndex(Index) /),   & 
     430                               count=(/ncell,size(field(1)%rval4d,2),1 /)) 
     431            DEALLOCATE(field_val3d) 
     432          ENDDO          
     433        ENDIF 
     434         
     435!        PRINT *,NF90_STRERROR(status) 
     436!        ncell=ncell+n 
     437 
     438!     ENDDO 
     439      
     440     ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
     441 
     442        ncell=0 
     443        DO ind=ind_b,ind_e 
     444          d=>domain_type(ind) 
     445   
     446          DO j=d%jj_begin+1,d%jj_end 
     447            DO i=d%ii_begin,d%ii_end-1 
     448              ncell=ncell+1 
     449            ENDDO 
     450          ENDDO 
     451 
     452          DO j=d%jj_begin,d%jj_end-1 
     453            DO i=d%ii_begin+1,d%ii_end 
     454              ncell=ncell+1 
     455            ENDDO 
     456          ENDDO 
     457        ENDDO 
     458 
     459        IF (field(1)%ndim==2) THEN 
     460          ALLOCATE(Field_val2d(ncell)) 
     461 
     462          n=0 
     463           
     464          DO ind=ind_b,ind_e 
     465            d=>domain_type(ind) 
     466            DO j=d%jj_begin+1,d%jj_end 
     467              DO i=d%ii_begin,d%ii_end-1 
     468                n=n+1 
     469                k=d%iim*(j-1)+i 
     470                Field_val2d(n)=field(ind)%rval2d(k+d%z_down) 
     471              ENDDO 
     472            ENDDO 
     473 
     474            DO j=d%jj_begin,d%jj_end-1 
     475              DO i=d%ii_begin+1,d%ii_end 
     476                n=n+1 
     477                k=d%iim*(j-1)+i 
     478                Field_val2d(n)=field(ind)%rval2d(k+d%z_up) 
     479              ENDDO 
     480            ENDDO 
     481          ENDDO 
     482           
     483          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       & 
     484                              Field_val2d,start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /)) 
     485          DEALLOCATE(field_val2d) 
     486 
     487        ELSE IF (field(1)%ndim==3) THEN 
     488          ALLOCATE(Field_val3d(ncell,size(field(1)%rval3d,2))) 
     489          n=0 
     490          DO ind=ind_b,ind_e 
     491            d=>domain_type(ind) 
     492            DO j=d%jj_begin+1,d%jj_end 
     493              DO i=d%ii_begin,d%ii_end-1 
     494                n=n+1 
     495                k=d%iim*(j-1)+i 
     496                Field_val3d(n,:)=field(ind)%rval3d(k+d%z_down,:) 
     497              ENDDO 
     498            ENDDO 
     499 
     500            DO j=d%jj_begin,d%jj_end-1 
     501              DO i=d%ii_begin+1,d%ii_end 
     502                n=n+1 
     503                k=d%iim*(j-1)+i 
     504                Field_val3d(n,:)=field(ind)%rval3d(k+d%z_up,:) 
     505              ENDDO 
     506            ENDDO 
     507          ENDDO 
     508           
     509          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1,FieldIndex(Index) /),   & 
     510                              count=(/ncell,size(field(1)%rval3d,2),1 /)) 
     511          DEALLOCATE(field_val3d) 
     512         
     513        ELSE IF (field(1)%ndim==4) THEN 
     514 
     515          DO q=1,FieldVarId(index)%size 
     516            ALLOCATE(Field_val3d(ncell,size(field(1)%rval4d,2))) 
     517            n=0 
     518            DO ind=ind_b,ind_e 
     519              d=>domain_type(ind) 
     520              DO j=d%jj_begin+1,d%jj_end 
     521                DO i=d%ii_begin,d%ii_end-1 
     522                  n=n+1 
     523                  k=d%iim*(j-1)+i 
     524                  Field_val3d(n,:)=field(ind)%rval4d(k+d%z_down,:,q) 
     525                ENDDO 
     526              ENDDO 
     527 
     528              DO j=d%jj_begin,d%jj_end-1 
     529                DO i=d%ii_begin+1,d%ii_end 
     530                  n=n+1 
     531                  k=d%iim*(j-1)+i 
     532                  Field_val3d(n,:)=field(ind)%rval4d(k+d%z_up,:,q) 
     533                ENDDO 
     534              ENDDO 
     535            ENDDO 
     536 
     537            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,1,FieldIndex(Index) /),   & 
     538                                count=(/ncell,size(field(1)%rval4d,2),1 /)) 
     539            DEALLOCATE(field_val3d) 
     540          ENDDO 
     541        ENDIF 
     542         
     543!        PRINT *,NF90_STRERROR(status) 
     544!        ncell=ncell+n 
     545! 
     546!     ENDDO 
     547      
     548     ENDIF 
     549     status=NF90_SYNC(FieldId(Index)) 
     550       
     551    END SUBROUTINE Writefield_gen 
     552 
     553       
     554 
     555    SUBROUTINE Writefield_mpi(name_in,field,nind) 
     556    USE netcdf_mod 
    81557    USE domain_mod 
    82558    use field_mod 
    … …  
    92568      TYPE(t_domain),POINTER :: d 
    93569      INTEGER :: Index 
    94       INTEGER :: ind,i,j,k,n,ncell,q 
     570      INTEGER :: ind,i,j,l,k,n,ncell,q 
    95571      INTEGER :: iie,jje,iin,jjn 
    96572      INTEGER :: status 
    … …  
    100576      INTEGER :: halo_size 
    101577      LOGICAL :: single 
     578      INTEGER :: displ 
    102579       
    103580       
    … …  
    120597      Index=GetFieldIndex(name) 
    121598      if (Index==-1) then 
    122         call create_header(name,field,nind) 
     599        call create_header_mpi(name,field,nind) 
    123600        Index=GetFieldIndex(name) 
    124601      else 
    … …  
    142619        ENDDO 
    143620 
     621        displ=FieldVarId(index)%displ 
     622         
    144623        IF (field(ind)%ndim==2) THEN 
    145624          ALLOCATE(Field_val2d(n)) 
    … …  
    155634          ENDDO 
    156635          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  & 
    157                               start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
     636                              start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /)) 
    158637          DEALLOCATE(field_val2d) 
    159638        ELSE IF (field(ind)%ndim==3) THEN 
    … …  
    169648            ENDDO 
    170649          ENDDO 
    171            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
    172                                count=(/n,size(field(1)%rval3d,2),1 /)) 
     650!           DO l=1,size(field(ind)%rval3d,2) 
     651            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ displ+ncell,1,FieldIndex(Index) /),   & 
     652                               count=(/n,size(field(ind)%rval3d,2),1 /)) 
     653!           ENDDO 
    173654          DEALLOCATE(field_val3d) 
    174655        ELSE IF (field(1)%ndim==4) THEN 
    … …  
    186667                ENDIF 
    187668              ENDDO 
    188             ENDDO 
    189             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
    190                                count=(/n,size(field(1)%rval4d,2),1 /)) 
     669             ENDDO 
     670!            DO l=1,size(field(ind)%rval4d,2) 
     671            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d(:,l),start=(/ displ+ncell,l,FieldIndex(Index) /),   & 
     672                               count=(/n,size(field(ind)%rval4d,2),1 /)) 
     673!            ENDDO 
    191674            DEALLOCATE(field_val3d) 
    192675          ENDDO          
    … …  
    219702          ENDDO 
    220703 
     704        displ=FieldVarId(index)%displ 
     705 
    221706        IF (field(ind)%ndim==2) THEN 
    222707          ALLOCATE(Field_val2d(n)) 
    … …  
    240725 
    241726          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       & 
    242                               Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
     727                              Field_val2d,start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /)) 
    243728          DEALLOCATE(field_val2d) 
    244729 
    … …  
    261746            ENDDO 
    262747          ENDDO 
    263            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
    264                                count=(/n,size(field(1)%rval3d,2),1 /)) 
     748!           DO l=1,size(field(ind)%rval3d,2) 
     749           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ displ+ncell,1,FieldIndex(Index) /),   & 
     750                               count=(/n,size(field(ind)%rval3d,2),1 /)) 
     751!           ENDDO 
    265752          DEALLOCATE(field_val3d) 
    266753        ELSE IF (field(1)%ndim==4) THEN 
    … …  
    285772            ENDDO 
    286773 
    287             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,1,FieldIndex(Index) /),   & 
    288                                 count=(/n,size(field(1)%rval4d,2),1 /)) 
     774!            DO l=1,size(field(ind)%rval4d,2) 
     775 
     776            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ displ+ncell,1,FieldIndex(Index) /),   & 
     777                                count=(/n,size(field(ind)%rval4d,2),1 /)) 
     778!            ENDDO 
    289779            DEALLOCATE(field_val3d) 
    290780          ENDDO 
    … …  
    299789     status=NF90_SYNC(FieldId(Index)) 
    300790       
    301     END SUBROUTINE Writefield 
    302        
     791    END SUBROUTINE Writefield_mpi 
     792     
     793     
    303794            
    304     SUBROUTINE Create_header(name,field,nind) 
    305     USE netcdf 
     795!   SUBROUTINE Create_header(name,field,nind) 
     796!   USE netcdf 
     797!   USE field_mod 
     798!   USE domain_mod 
     799!   USE spherical_geom_mod 
     800!   USE dimensions 
     801!   USE geometry 
     802!   IMPLICIT NONE 
     803!     CHARACTER(LEN=*) :: name 
     804!     TYPE(t_field),POINTER :: field(:) 
     805!     INTEGER,OPTIONAL,INTENT(IN) :: nind 
     806!     INTEGER :: ncell 
     807!     INTEGER :: nvert 
     808!     REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:) 
     809!     TYPE(t_domain),POINTER :: d 
     810!     INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId 
     811!     INTEGER :: dim3id,dim4id 
     812!     INTEGER :: status 
     813!     INTEGER :: ind,i,j,k,n,q 
     814!     INTEGER :: iie,jje,iin,jjn 
     815!     INTEGER :: ind_b,ind_e 
     816!     INTEGER :: halo_size 
     817!     LOGICAL :: single  
     818!     INTEGER :: nij 
     819!          
     820!     NbField=NbField+1 
     821!     FieldName(NbField)=TRIM(ADJUSTL(name)) 
     822!     FieldIndex(NbField)=1 
     823!      
     824!     IF (PRESENT(nind)) THEN 
     825!       ind_b=nind 
     826!       ind_e=nind 
     827!       halo_size=1 
     828!       single=.TRUE. 
     829!     ELSE 
     830!       ind_b=1 
     831!       ind_e=ndomain 
     832!       halo_size=0 
     833!       single=.FALSE. 
     834!     ENDIF 
     835!      
     836!     ncell=0 
     837!      
     838!     IF (Field(ind_b)%field_type==field_T) THEN 
     839!       nvert=6 
     840!        
     841!       DO ind=ind_b,ind_e 
     842!         d=>domain(ind) 
     843!        
     844!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     845!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     846!             IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1 
     847!           ENDDO 
     848!         ENDDO 
     849 
     850!       END DO 
     851!      
     852!       status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
     853!       FieldId(NbField)=ncid 
     854!       status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
     855!       status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
     856 
     857!       IF (Field(ind_b)%ndim==2)  THEN 
     858!         FieldVarId(NbField)%size=1 
     859!         ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     860!       ELSE IF (Field(ind_b)%ndim==3)  THEN 
     861!         FieldVarId(NbField)%size=1 
     862!         ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     863!         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) 
     864!       ELSE IF (Field(1)%ndim==4) THEN 
     865!         FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
     866!         ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
     867!         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) 
     868!          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id) 
     869!       ENDIF 
     870!      
     871!       status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
     872!      
     873!       status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) 
     874!       status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
     875!       status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
     876!       status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") 
     877!       status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) 
     878!       status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
     879!       status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
     880!       status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") 
     881!       status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
     882!       status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
     883 
     884!       IF (Field(ind_b)%ndim==2) THEN 
     885!         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
     886!         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     887!       ELSE IF (Field(ind_b)%ndim==3) THEN 
     888!         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
     889!         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     890!       ELSE IF (Field(ind_b)%ndim==4) THEN 
     891!         DO i=1,FieldVarId(NbField)%size 
     892!           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),  & 
     893!                                 FieldVarId(NbField)%nc_id(i)) 
     894!           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat") 
     895!         ENDDO         
     896!       ENDIF  
     897!  
     898!      
     899!       status = NF90_ENDDEF(ncid)       
     900 
     901!       ncell=1 
     902!       DO ind=ind_b,ind_e 
     903!         d=>domain(ind) 
     904!  
     905!         n=0 
     906!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     907!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     908!             IF (single .OR. d%own(i,j)) n=n+1 
     909!           ENDDO 
     910!         ENDDO 
     911!          
     912!        ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n)) 
     913!          
     914!         n=0   
     915!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     916!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     917!               IF (d%own(i,j) .OR. single) THEN 
     918!               n=n+1 
     919!               CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n)) 
     920!               lon(n)=lon(n)*180/Pi 
     921!               lat(n)=lat(n)*180/Pi 
     922!               DO k=0,5 
     923!                 CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n)) 
     924!                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
     925!                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     926!               ENDDO 
     927!             ENDIF 
     928!           ENDDO 
     929!         ENDDO 
     930!         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /)) 
     931!         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /)) 
     932!         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
     933!         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
     934!  
     935!         ncell=ncell+n 
     936!         DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
     937!     END DO  
     938 
     939!   ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
     940!       nvert=3 
     941!       DO ind=ind_b,ind_e 
     942!         d=>domain(ind) 
     943!        
     944!         DO j=d%jj_begin+1,d%jj_end 
     945!           DO i=d%ii_begin,d%ii_end-1 
     946!             ncell=ncell+1 
     947!           ENDDO 
     948!         ENDDO 
     949 
     950!         DO j=d%jj_begin,d%jj_end-1 
     951!           DO i=d%ii_begin+1,d%ii_end 
     952!             ncell=ncell+1 
     953!           ENDDO 
     954!         ENDDO 
     955 
     956!       END DO 
     957!      
     958!       status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
     959!       FieldId(NbField)=ncid 
     960!       status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
     961!       status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
     962 
     963!       IF (Field(ind_b)%ndim==2)  THEN 
     964!         FieldVarId(NbField)%size=1 
     965!         ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     966!       ELSE IF (Field(ind_b)%ndim==3)  THEN 
     967!         FieldVarId(NbField)%size=1 
     968!         ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     969!         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) 
     970!       ELSE IF (Field(1)%ndim==4) THEN 
     971!         FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
     972!         ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
     973!         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) 
     974!       ENDIF 
     975 
     976 
     977!      
     978!       status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
     979!      
     980!       status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) 
     981!       status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
     982!       status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
     983!       status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") 
     984!       status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) 
     985!       status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
     986!       status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
     987!       status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") 
     988!       status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
     989!       status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
     990 
     991 
     992!       IF (Field(ind_b)%ndim==2) THEN 
     993!         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
     994!         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     995!       ELSE IF (Field(ind_b)%ndim==3) THEN 
     996!         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
     997!         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     998!       ELSE IF (Field(ind_b)%ndim==4) THEN 
     999!         DO q=1,FieldVarId(NbField)%size 
     1000!           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),NF90_DOUBLE,             & 
     1001!                                 (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q)) 
     1002!           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat") 
     1003!         ENDDO         
     1004!       ENDIF  
     1005!        
     1006!       status = NF90_ENDDEF(ncid)       
     1007 
     1008!       ncell=1 
     1009!       DO ind=ind_b,ind_e 
     1010!         d=>domain(ind) 
     1011!         CALL swap_geometry(ind) 
     1012!         CALL swap_dimensions(ind) 
     1013!  
     1014!         n=0 
     1015!         DO j=jj_begin+1,jj_end 
     1016!           DO i=ii_begin,ii_end-1 
     1017!             n=n+1 
     1018!           ENDDO 
     1019!         ENDDO 
     1020 
     1021!         DO j=jj_begin,jj_end-1 
     1022!           DO i=ii_begin+1,ii_end 
     1023!             n=n+1 
     1024!           ENDDO 
     1025!         ENDDO 
     1026 
     1027!        ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n)) 
     1028!          
     1029!         n=0   
     1030!       
     1031!         DO j=jj_begin+1,jj_end 
     1032!           DO i=ii_begin,ii_end-1 
     1033!             nij=(j-1)*iim+i 
     1034!             n=n+1 
     1035!             CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n)) 
     1036!             lon(n)=lon(n)*180/Pi 
     1037!!             IF (lon(n)<0) lon(n)=lon(n)+360 
     1038!             lat(n)=lat(n)*180/Pi 
     1039!             CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
     1040!             CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) 
     1041!             CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) 
     1042 
     1043!             DO k=0,2 
     1044!               bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
     1045!               bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     1046!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 
     1047!             ENDDO 
     1048!           ENDDO 
     1049!         ENDDO 
     1050 
     1051!         DO j=jj_begin,jj_end-1 
     1052!           DO i=ii_begin+1,ii_end 
     1053!             nij=(j-1)*iim+i 
     1054!             n=n+1 
     1055!             CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n)) 
     1056!             lon(n)=lon(n)*180/Pi 
     1057!              IF (lon(n)<0) lon(n)=lon(n)+360 
     1058!             lat(n)=lat(n)*180/Pi 
     1059!             CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
     1060!             CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) 
     1061!             CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) 
     1062 
     1063!             DO k=0,2 
     1064!               bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
     1065!               bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     1066!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 
     1067!             ENDDO 
     1068!           ENDDO 
     1069!         ENDDO 
     1070!          
     1071!          
     1072!         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /)) 
     1073!         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /)) 
     1074!         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
     1075!         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
     1076!  
     1077!         ncell=ncell+n 
     1078!         DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
     1079!     END DO           
     1080!   ENDIF 
     1081 
     1082 
     1083!      
     1084!      status = NF90_CLOSE(ncid) 
     1085 
     1086!   END SUBROUTINE Create_Header 
     1087 
     1088 
     1089 
     1090    
     1091    SUBROUTINE Create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in) 
     1092    USE netcdf_mod 
     1093    USE field_mod 
     1094    USE domain_mod 
     1095    USE metric 
     1096    USE spherical_geom_mod 
     1097!    USE dimensions 
     1098!    USE geometry 
     1099    IMPLICIT NONE 
     1100      CHARACTER(LEN=*),INTENT(IN) :: name_in 
     1101      TYPE(t_field),POINTER :: field(:) 
     1102      TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:) 
     1103      INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in 
     1104      INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in 
     1105      INTEGER :: ncell 
     1106      INTEGER :: nvert 
     1107      REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:) 
     1108      TYPE(t_domain),POINTER :: d 
     1109      INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId 
     1110      INTEGER :: dim3id,dim4id 
     1111      INTEGER :: status 
     1112      INTEGER :: ind,i,j,k,n,q 
     1113      INTEGER :: iie,jje,iin,jjn 
     1114      INTEGER :: ind_b,ind_e 
     1115      INTEGER :: halo_size 
     1116      LOGICAL :: single  
     1117      INTEGER :: nij 
     1118      CHARACTER(LEN=255) :: name 
     1119       
     1120             
     1121      name=TRIM(ADJUSTL(name_in)) 
     1122 
     1123      IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN 
     1124        name=TRIM(name)//"_"//TRIM(int2str(ind_b)) 
     1125        ind_b=ind_b_in 
     1126        ind_e=ind_b_in 
     1127        halo_size=1 
     1128        single=.TRUE. 
     1129      ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
     1130        name=TRIM(name)//"_"//TRIM(int2str(ind_e)) 
     1131        ind_b=ind_e_in 
     1132        ind_e=ind_e_in 
     1133        halo_size=1 
     1134        single=.TRUE. 
     1135      ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
     1136        ind_b=ind_b_in 
     1137        ind_e=ind_e_in 
     1138        halo_size=0 
     1139        single=.FALSE. 
     1140      ELSE 
     1141        halo_size=0 
     1142        ind_b=1 
     1143        ind_e=ndomain 
     1144        single=.FALSE. 
     1145      ENDIF       
     1146                 
     1147      NbField=NbField+1 
     1148      FieldName(NbField)=TRIM(ADJUSTL(name)) 
     1149      FieldIndex(NbField)=1 
     1150       
     1151      ncell=0 
     1152       
     1153      IF (Field(ind_b)%field_type==field_T) THEN 
     1154        nvert=6 
     1155         
     1156        DO ind=ind_b,ind_e 
     1157          d=>domain_type(ind) 
     1158         
     1159          DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     1160            DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     1161              IF (single .OR. d%assign_domain(i,j)==ind) ncell=ncell+1 
     1162            ENDDO 
     1163          ENDDO 
     1164 
     1165        END DO 
     1166       
     1167        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
     1168        FieldId(NbField)=ncid 
     1169        status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
     1170        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
     1171 
     1172        IF (Field(ind_b)%ndim==2)  THEN 
     1173          FieldVarId(NbField)%size=1 
     1174          ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     1175        ELSE IF (Field(ind_b)%ndim==3)  THEN 
     1176          FieldVarId(NbField)%size=1 
     1177          ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     1178          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) 
     1179        ELSE IF (Field(1)%ndim==4) THEN 
     1180          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
     1181          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
     1182          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) 
     1183!          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id) 
     1184        ENDIF 
     1185       
     1186        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
     1187       
     1188        status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) 
     1189        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
     1190        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
     1191        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") 
     1192        status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) 
     1193        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
     1194        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
     1195        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") 
     1196        status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
     1197        status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
     1198 
     1199        IF (Field(ind_b)%ndim==2) THEN 
     1200          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
     1201          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1202        ELSE IF (Field(ind_b)%ndim==3) THEN 
     1203          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
     1204          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1205        ELSE IF (Field(ind_b)%ndim==4) THEN 
     1206          DO i=1,FieldVarId(NbField)%size 
     1207            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),  & 
     1208                                  FieldVarId(NbField)%nc_id(i)) 
     1209            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat") 
     1210          ENDDO         
     1211        ENDIF  
     1212   
     1213       
     1214        status = NF90_ENDDEF(ncid)       
     1215 
     1216!        ncell=1 
     1217!        DO ind=ind_b,ind_e 
     1218!          d=>domain_type(ind) 
     1219   
     1220!          n=0 
     1221!          DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     1222!            DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     1223!              IF (single .OR. d%assign_domain(i,j)==ind) n=n+1 
     1224!            ENDDO 
     1225!          ENDDO 
     1226           
     1227         ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell)) 
     1228           
     1229         n=0   
     1230         DO ind=ind_b,ind_e 
     1231           d=>domain_type(ind) 
     1232           DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     1233             DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     1234               IF (d%assign_domain(i,j)==ind .OR. single) THEN 
     1235                 n=n+1 
     1236                 CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n)) 
     1237                 lon(n)=lon(n)*180/Pi 
     1238                 lat(n)=lat(n)*180/Pi 
     1239                 DO k=0,5 
     1240                   CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n)) 
     1241                   bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
     1242                   bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     1243                 ENDDO 
     1244               ENDIF 
     1245             ENDDO 
     1246           ENDDO 
     1247         ENDDO 
     1248         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /)) 
     1249         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /)) 
     1250         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
     1251         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
     1252   
     1253!         ncell=ncell+n 
     1254         DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
     1255!      END DO  
     1256 
     1257    ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
     1258        nvert=3 
     1259        DO ind=ind_b,ind_e 
     1260          d=>domain_type(ind) 
     1261         
     1262          DO j=d%jj_begin+1,d%jj_end 
     1263            DO i=d%ii_begin,d%ii_end-1 
     1264              ncell=ncell+1 
     1265            ENDDO 
     1266          ENDDO 
     1267 
     1268          DO j=d%jj_begin,d%jj_end-1 
     1269            DO i=d%ii_begin+1,d%ii_end 
     1270              ncell=ncell+1 
     1271            ENDDO 
     1272          ENDDO 
     1273 
     1274        END DO 
     1275       
     1276        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
     1277        FieldId(NbField)=ncid 
     1278        status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
     1279        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
     1280 
     1281        IF (Field(ind_b)%ndim==2)  THEN 
     1282          FieldVarId(NbField)%size=1 
     1283          ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     1284        ELSE IF (Field(ind_b)%ndim==3)  THEN 
     1285          FieldVarId(NbField)%size=1 
     1286          ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     1287          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) 
     1288        ELSE IF (Field(1)%ndim==4) THEN 
     1289          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
     1290          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
     1291          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) 
     1292        ENDIF 
     1293 
     1294 
     1295       
     1296        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
     1297       
     1298        status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) 
     1299        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
     1300        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
     1301        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") 
     1302        status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) 
     1303        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
     1304        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
     1305        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") 
     1306        status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
     1307        status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
     1308 
     1309 
     1310        IF (Field(ind_b)%ndim==2) THEN 
     1311          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
     1312          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1313        ELSE IF (Field(ind_b)%ndim==3) THEN 
     1314          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
     1315          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1316        ELSE IF (Field(ind_b)%ndim==4) THEN 
     1317          DO q=1,FieldVarId(NbField)%size 
     1318            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),NF90_DOUBLE,             & 
     1319                                  (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q)) 
     1320            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat") 
     1321          ENDDO         
     1322        ENDIF  
     1323         
     1324        status = NF90_ENDDEF(ncid)       
     1325 
     1326!        ncell=1 
     1327!        DO ind=ind_b,ind_e 
     1328!          d=>domain_type(ind) 
     1329!   
     1330!          n=0 
     1331!          DO j=d%jj_begin+1,d%jj_end 
     1332!            DO i=d%ii_begin,d%ii_end-1 
     1333!              n=n+1 
     1334!            ENDDO 
     1335!          ENDDO 
     1336! 
     1337!          DO j=d%jj_begin,d%jj_end-1 
     1338!            DO i=d%ii_begin+1,d%ii_end 
     1339!              n=n+1 
     1340!            ENDDO 
     1341!          ENDDO 
     1342 
     1343!         ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n)) 
     1344         ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell)) 
     1345           
     1346         n=0   
     1347        
     1348         DO ind=ind_b,ind_e 
     1349           d=>domain_type(ind) 
     1350           DO j=d%jj_begin+1,d%jj_end 
     1351             DO i=d%ii_begin,d%ii_end-1 
     1352               nij=(j-1)*d%iim+i 
     1353               n=n+1 
     1354               CALL xyz2lonlat(d%vertex(:,vdown,i,j),lon(n),lat(n)) 
     1355               lon(n)=lon(n)*180/Pi 
     1356 !              IF (lon(n)<0) lon(n)=lon(n)+360 
     1357               lat(n)=lat(n)*180/Pi 
     1358!               CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
     1359!               CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) 
     1360!               CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) 
     1361 
     1362               CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n)) 
     1363               CALL xyz2lonlat(d%xyz(:,i,j-1),bounds_lon(1,n), bounds_lat(1,n)) 
     1364               CALL xyz2lonlat(d%xyz(:,i+1,j-1),bounds_lon(2,n), bounds_lat(2,n)) 
     1365 
     1366               DO k=0,2 
     1367                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
     1368                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     1369!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 
     1370               ENDDO 
     1371             ENDDO 
     1372           ENDDO 
     1373 
     1374           DO j=d%jj_begin,d%jj_end-1 
     1375             DO i=d%ii_begin+1,d%ii_end 
     1376               nij=(j-1)*d%iim+i 
     1377               n=n+1 
     1378!              CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n)) 
     1379               CALL xyz2lonlat(d%vertex(:,vup,i,j),lon(n),lat(n)) 
     1380               lon(n)=lon(n)*180/Pi 
     1381!              IF (lon(n)<0) lon(n)=lon(n)+360 
     1382               lat(n)=lat(n)*180/Pi 
     1383!              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
     1384!              CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) 
     1385!              CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) 
     1386               CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n)) 
     1387               CALL xyz2lonlat(d%xyz(:,i,j+1),bounds_lon(1,n), bounds_lat(1,n)) 
     1388               CALL xyz2lonlat(d%xyz(:,i-1,j+1),bounds_lon(2,n), bounds_lat(2,n)) 
     1389 
     1390               DO k=0,2 
     1391                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
     1392                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     1393!                 IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 
     1394               ENDDO 
     1395             ENDDO 
     1396           ENDDO 
     1397         ENDDO  
     1398           
     1399         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /)) 
     1400         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /)) 
     1401         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
     1402         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
     1403   
     1404!          ncell=ncell+n 
     1405          DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
     1406!      END DO           
     1407    ENDIF 
     1408 
     1409 
     1410       
     1411!      status = NF90_CLOSE(ncid) 
     1412 
     1413    END SUBROUTINE Create_Header_gen 
     1414 
     1415    SUBROUTINE Create_header_mpi(name,field,nind) 
     1416    USE netcdf_mod 
    3061417    USE field_mod 
    3071418    USE domain_mod 
    … …  
    3091420    USE dimensions 
    3101421    USE geometry 
     1422    USE mpi_mod 
     1423    USE mpipara 
    3111424    IMPLICIT NONE 
    3121425      CHARACTER(LEN=*) :: name 
    … …  
    3261439      LOGICAL :: single  
    3271440      INTEGER :: nij 
     1441      INTEGER :: ncell_glo(0:mpi_size-1) 
     1442      INTEGER :: displ, ncell_tot 
     1443       
    3281444           
    3291445      NbField=NbField+1 
    … …  
    3591475        END DO 
    3601476       
    361         status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
     1477        CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr) 
     1478 
     1479        displ=0 
     1480        DO i=1,mpi_rank 
     1481          displ=displ+ncell_glo(i-1) 
     1482        ENDDO 
     1483        FieldVarId(NbField)%displ=displ 
     1484        ncell_tot=sum(ncell_glo(:)) 
     1485         
     1486        status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc', IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid) 
    3621487        FieldId(NbField)=ncid 
    363         status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
     1488        status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId) 
    3641489        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
    3651490 
    … …  
    3941519          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
    3951520          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1521          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/)) 
    3961522        ELSE IF (Field(ind_b)%ndim==3) THEN 
    3971523          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
    3981524          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1525          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,size(field(ind_b)%rval3d,2),1/)) 
    3991526        ELSE IF (Field(ind_b)%ndim==4) THEN 
    4001527          DO i=1,FieldVarId(NbField)%size 
    … …  
    4021529                                  FieldVarId(NbField)%nc_id(i)) 
    4031530            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat") 
     1531            status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED, (/ncell_tot,size(field(ind_b)%rval4d,2),1/)) 
    4041532          ENDDO         
    4051533        ENDIF  
    … …  
    4371565            ENDDO 
    4381566          ENDDO 
    439           status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /)) 
    440           status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /)) 
    441           status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
    442           status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
     1567          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /)) 
     1568          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /)) 
     1569          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) 
     1570          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) 
    4431571   
    4441572          ncell=ncell+n 
    … …  
    4641592 
    4651593        END DO 
    466        
    467         status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
     1594         
     1595        CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr) 
     1596 
     1597        displ=0 
     1598        DO i=1,mpi_rank 
     1599          displ=displ+ncell_glo(i-1) 
     1600        ENDDO 
     1601        FieldVarId(NbField)%displ=displ 
     1602        ncell_tot=sum(ncell_glo(:)) 
     1603               
     1604        status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc',IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid) 
     1605!        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
    4681606        FieldId(NbField)=ncid 
    469         status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
     1607        status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId) 
    4701608        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
    4711609 
    … …  
    5021640          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
    5031641          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1642          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/)) 
    5041643        ELSE IF (Field(ind_b)%ndim==3) THEN 
    5051644          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
    5061645          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1646          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,size(field(ind_b)%rval3d,2),1/)) 
    5071647        ELSE IF (Field(ind_b)%ndim==4) THEN 
    5081648          DO q=1,FieldVarId(NbField)%size 
    … …  
    5101650                                  (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q)) 
    5111651            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat") 
     1652           status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED, (/ncell_tot,size(field(ind_b)%rval4d,2),1/)) 
    5121653          ENDDO         
    5131654        ENDIF  
    … …  
    5791720           
    5801721           
    581           status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /)) 
    582           status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /)) 
    583           status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
    584           status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
     1722          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /)) 
     1723          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /)) 
     1724          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) 
     1725          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) 
    5851726   
    5861727          ncell=ncell+n 
    … …  
    5931734!      status = NF90_CLOSE(ncid) 
    5941735 
    595     end subroutine Create_Header 
     1736    end subroutine Create_Header_mpi   
    5961737    
    597    
     1738   SUBROUTINE Close_files 
     1739   USE netcdf 
     1740   IMPLICIT NONE 
     1741     INTEGER :: i,k,status 
     1742      
     1743     DO i=1,NbField 
     1744         status=NF90_CLOSE(FieldId(i)) 
     1745     ENDDO 
     1746   END SUBROUTINE  Close_files 
     1747      
     1748      
    5981749  function int2str(int) 
    5991750    implicit none 
Note: See TracChangeset for help on using the changeset viewer.