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

Implementation of parallelism
implementation of iterative laplacian for dissipation

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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 
Note: See TracChangeset for help on using the changeset viewer.