Changeset 151

Timestamp:

05/13/13 14:30:31 (11 years ago)

Author:

ymipsl

Message:

Implementation of mixte parallelism MPI/OpenMP into src directory

YM

Location:

codes/icosagcm/trunk/src

Files:

• advect.f90 (modified) (10 diffs)
• advect_tracer.f90 (modified) (12 diffs)
• caldyn.f90 (modified) (4 diffs)
• caldyn_gcm.f90 (modified) (31 diffs)
• caldyn_gcm_opt.f90 (deleted)
• check_conserve.f90 (modified) (9 diffs)
• check_conserve.f90~ (deleted)
• diffff (deleted)
• dimensions.f90 (modified) (2 diffs)
• dissip_gcm.f90 (modified) (35 diffs)
• domain.f90 (modified) (6 diffs)
• dynetat0_gcm_mod.f90 (modified) (2 diffs)
• dynetat0_hz_mod.f90 (modified) (3 diffs)
• exner.f90 (modified) (6 diffs)
• geometry.f90 (modified) (8 diffs)
• icosa_gcm.f90 (modified) (3 diffs)
• mpi_mod.F90 (modified) (3 diffs)
• mpipara.F90 (modified) (2 diffs)
• omp_para.F90 (added)
• pression.f90 (modified) (1 diff)
• theta_rhodz.f90 (modified) (5 diffs)
• time.f90 (modified) (2 diffs)
• timeloop_gcm.f90 (modified) (23 diffs)
• trace.F90 (modified) (3 diffs)
• transfert.F90 (modified) (1 diff)
• transfert_mpi.f90 (modified) (15 diffs)
• wind.f90 (modified) (2 diffs)
• write_field.f90 (modified) (23 diffs)

codes/icosagcm/trunk/src/advect.f90

@@ -53 +53 @@
   SUBROUTINE compute_gradq3d(qi,one_over_sqrt_leng,gradq3d)
     USE trace
+    USE omp_para
     IMPLICIT NONE
     REAL(rstd),INTENT(IN)  :: qi(iim*jjm,llm)
@@ -72 +73 @@
     ! Compute gradient at triangles solving a linear system
     ! arr = area of triangle joining centroids of hexagons
-    Do l = 1,llm 
+     DO l = ll_begin,ll_end 
        DO j=jj_begin-1,jj_end+1
           DO i=ii_begin-1,ii_end+1
@@ -90 +91 @@
       
     DO k=1,3
-      DO l =1,llm
+      DO l =ll_begin,ll_end
         DO j=jj_begin,jj_end
            DO i=ii_begin,ii_end
@@ -103 +104 @@
 
     !============================================================================================= LIMITING 
-    ! GO TO 120 
-    DO l =1,llm
+    DO l =ll_begin,ll_end
        DO j=jj_begin,jj_end
           DO i=ii_begin,ii_end
@@ -110 +110 @@
              maggrd =  dot_product(gradq3d(n,l,:),gradq3d(n,l,:))
              maggrd = sqrt(maggrd) 
-!             leng = max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), &
-!                  sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2),  &
-!                  sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2))
-!             maxq_c = qi(n,l) + maggrd*sqrt(leng(n))
-!             minq_c = qi(n,l) - maggrd*sqrt(leng(n))
              maxq_c = qi(n,l) + maggrd*one_over_sqrt_leng(n)
              minq_c = qi(n,l) - maggrd*one_over_sqrt_leng(n)
@@ -137 +132 @@
   SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc)
     USE trace
+    USE omp_para
     IMPLICIT NONE
     REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3)
@@ -152 +148 @@
     
     ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau
-    DO l = 1,llm
+    DO l = ll_begin,ll_end
        DO j=jj_begin-1,jj_end+1
           DO i=ii_begin-1,ii_end+1
@@ -213 +209 @@
   SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi) 
     USE trace
+    USE omp_para
     IMPLICIT NONE
     LOGICAL, INTENT(IN)       :: update_mass
@@ -230 +227 @@
     ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon
     ! ne*hfluxt>0 iff outgoing
-    DO l = 1,llm
+    DO l = ll_begin,ll_end
        DO j=jj_begin-1,jj_end+1
           DO i=ii_begin-1,ii_end+1
@@ -265 +262 @@
 
     ! update q and, if update_mass, update mass
-    DO l = 1,llm
+    DO l = ll_begin,ll_end
        DO j=jj_begin,jj_end
           DO i=ii_begin,ii_end

codes/icosagcm/trunk/src/advect_tracer.f90

@@ -9 +9 @@
   TYPE(t_field),POINTER :: f_cc(:)  ! starting point of backward-trajectory (Miura approach)
   TYPE(t_field),POINTER :: f_one_over_sqrt_leng(:)
- 
+
+  TYPE(t_message) :: req_u, req_wfluxt, req_q, req_rhodz, req_gradq3d
+
   REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz
+
+! temporary shared variable for vlz
+  TYPE(t_field),POINTER :: f_dzqw(:)   ! vertical finite difference of q 
+  TYPE(t_field),POINTER :: f_adzqw(:)  ! abs(dzqw)
+  TYPE(t_field),POINTER :: f_dzq(:)    ! limited slope of q
+  TYPE(t_field),POINTER :: f_wq(:)     ! time-integrated flux of q
 
   PUBLIC init_advect_tracer, advect_tracer
@@ -28 +36 @@
     CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc')
     CALL allocate_field(f_one_over_sqrt_leng,field_t,type_real, name='one_over_sqrt_leng')
-
+    CALL allocate_field(f_dzqw, field_t, type_real, llm, name='dzqw')
+    CALL allocate_field(f_adzqw, field_t, type_real, llm, name='adzqw')
+    CALL allocate_field(f_dzq, field_t, type_real, llm, name='dzq')
+    CALL allocate_field(f_wq, field_t, type_real, llm+1, name='wq')
+    
     DO ind=1,ndomain
        CALL swap_dimensions(ind)
@@ -56 +68 @@
     REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:)
     REAL(rstd),POINTER :: rhodz(:,:), u(:,:) 
-    INTEGER :: ind,k
+! temporary shared variable for vlz
+    REAL(rstd),POINTER ::  dzqw(:,:)         ! vertical finite difference of q 
+    REAL(rstd),POINTER ::  adzqw(:,:)        ! abs(dzqw)
+    REAL(rstd),POINTER ::  dzq(:,:)          ! limited slope of q
+    REAL(rstd),POINTER ::  wq(:,:)           ! time-integrated flux of q
+    
+     INTEGER :: ind,k
+    LOGICAL,SAVE :: first=.TRUE.
+!$OMP THREADPRIVATE(first)
+
+    IF (first) THEN
+      first=.FALSE.
+      CALL init_message(f_u,req_e1_vect,req_u)
+      CALL init_message(f_wfluxt,req_i1,req_wfluxt)
+      CALL init_message(f_q,req_i1,req_q)
+      CALL init_message(f_rhodz,req_i1,req_rhodz)
+      CALL init_message(f_gradq3d,req_i1,req_gradq3d)
+    ENDIF
+    
+!$OMP BARRIER
 
     CALL trace_start("advect_tracer") 
 
-    CALL transfert_request(f_u,req_e1_vect)
-!    CALL transfert_request(f_hfluxt,req_e1)  ! BUG : This (unnecessary) transfer makes the computation go wrong
-    CALL transfert_request(f_wfluxt,req_i1)
-    CALL transfert_request(f_q,req_i1)
-    CALL transfert_request(f_rhodz,req_i1)
-        
-!    IF (is_mpi_root) PRINT *, 'Advection scheme'
-
-!    DO ind=1,ndomain
-!       CALL swap_dimensions(ind)
-!       CALL swap_geometry(ind)
-!       normal  = f_normal(ind)
-!       tangent = f_tangent(ind)
-!       cc      = f_cc(ind)
-!       u       = f_u(ind)
-!       q       = f_q(ind)
-!       rhodz   = f_rhodz(ind)
-!       hfluxt  = f_hfluxt(ind) 
-!       wfluxt  = f_wfluxt(ind) 
-!       gradq3d = f_gradq3d(ind)
-!
-!       ! 1/2 vertical transport
-!       DO k = 1, nqtot
-!          CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k))
-!       END DO
-!
-!       ! horizontal transport
-!       CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc) 
-!       DO k = 1,nqtot
-!          CALL compute_gradq3d(q(:,:,k),gradq3d)
-!          CALL compute_advect_horiz(k==nqtot,hfluxt,cc,gradq3d, rhodz,q(:,:,k))
-!       END DO
-!
-!       ! 1/2 vertical transport
-!       DO k = 1,nqtot
-!          CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k))
-!       END DO
-!    END DO
-
+    CALL send_message(f_u,req_u)
+    CALL send_message(f_wfluxt,req_wfluxt)
+    CALL send_message(f_q,req_q)
+    CALL send_message(f_rhodz,req_rhodz)
+    CALL wait_message(req_u)
+    CALL wait_message(req_wfluxt)
+    CALL wait_message(req_q)
+    CALL wait_message(req_rhodz)
+    
     ! 1/2 vertical transport + back-trajectories
     DO ind=1,ndomain
@@ -110 +111 @@
        rhodz   = f_rhodz(ind)
        wfluxt  = f_wfluxt(ind) 
+       dzqw    = f_dzqw(ind)
+       adzqw   = f_adzqw(ind)
+       dzq     = f_dzq(ind)
+       wq      = f_wq(ind) 
 
        DO k = 1, nqtot
-          CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k),1)
+          CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k),1,dzqw, adzqw, dzq, wq)
        END DO
 
        CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc) 
-    END DO
-
-!    CALL transfert_request(f_q,req_i1)      ! necessary ?
-!    CALL transfert_request(f_rhodz,req_i1)  ! necessary ?
+
+    END DO
+
 
     ! horizontal transport - split in two to place transfer of gradq3d
-
+!$OMP BARRIER
     DO k = 1, nqtot
        DO ind=1,ndomain
@@ -133 +137 @@
        END DO
 
-       CALL transfert_request(f_gradq3d,req_i1)
-
+       CALL send_message(f_gradq3d,req_gradq3d)
+       CALL wait_message(req_gradq3d)
 
 
@@ -149 +153 @@
     END DO 
     
-!    CALL transfert_request(f_q,req_i1)      ! necessary ?
-!    CALL transfert_request(f_rhodz,req_i1)  ! necessary ?
-    
     ! 1/2 vertical transport
+!$OMP BARRIER
+
     DO ind=1,ndomain
        CALL swap_dimensions(ind)
@@ -159 +162 @@
        rhodz   = f_rhodz(ind)
        wfluxt  = f_wfluxt(ind) 
+       dzqw    = f_dzqw(ind)
+       adzqw   = f_adzqw(ind)
+       dzq     = f_dzq(ind)
+       wq      = f_wq(ind) 
+
        DO k = 1,nqtot
-          CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k),0)
+          CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k),0, dzqw, adzqw, dzq, wq)
        END DO
+
     END DO
 
     CALL trace_end("advect_tracer")
 
+!$OMP BARRIER
+
   END SUBROUTINE advect_tracer
 
-  SUBROUTINE vlz(update_mass, fac,wfluxt,mass, q, halo)
+  SUBROUTINE vlz(update_mass, fac,wfluxt,mass, q, halo, dzqw, adzqw, dzq, wq)
     !
     !     Auteurs:   P.Le Van, F.Hourdin, F.Forget, T. Dubos
@@ -178 +189 @@
     !    ********************************************************************
     USE trace
+    USE omp_para
     IMPLICIT NONE
     LOGICAL, INTENT(IN)       :: update_mass
@@ -185 +197 @@
     INTEGER, INTENT(IN) :: halo
 
-    REAL(rstd) :: dq(iim*jjm,llm), & ! increase of q
-         dzqw(iim*jjm,llm),        & ! vertical finite difference of q 
-         adzqw(iim*jjm,llm),       & ! abs(dzqw)
-         dzq(iim*jjm,llm),         & ! limited slope of q
-         wq(iim*jjm,llm+1)           ! time-integrated flux of q
+! temporary shared variable
+    REAL(rstd),INTENT(INOUT) :: dzqw(iim*jjm,llm),        & ! vertical finite difference of q 
+                                adzqw(iim*jjm,llm),       & ! abs(dzqw)
+                                dzq(iim*jjm,llm),         & ! limited slope of q
+                                wq(iim*jjm,llm+1)           ! time-integrated flux of q
+
+
     REAL(rstd) :: dzqmax, newmass, sigw, qq, w
     INTEGER :: i,ij,l,j
@@ -196 +210 @@
 
     ! finite difference of q
-    DO l=2,llm
+
+     DO l=ll_beginp1,ll_end
        DO j=jj_begin-halo,jj_end+halo
           DO i=ii_begin-halo,ii_end+halo
@@ -206 +221 @@
     ENDDO
 
+!--> flush dzqw, adzqw
+!$OMP BARRIER
+
     ! minmod-limited slope of q
     ! dzq = slope*dz, i.e. the reconstructed q varies by dzq inside level l
-    DO l=2,llm-1
+
+     DO l=ll_beginp1,ll_endm1
        DO j=jj_begin-halo,jj_end+halo
           DO i=ii_begin-halo,ii_end+halo
@@ -223 +242 @@
     ENDDO
 
+
     ! 0 slope in top and bottom layers
-    DO j=jj_begin-halo,jj_end+halo
-       DO i=ii_begin-halo,ii_end+halo
-          ij=(j-1)*iim+i
-          dzq(ij,1)=0.
+    IF (omp_first) THEN
+      DO j=jj_begin-halo,jj_end+halo
+         DO i=ii_begin-halo,ii_end+halo
+           ij=(j-1)*iim+i
+           dzq(ij,1)=0.
+         ENDDO
+      ENDDO
+    ENDIF
+      
+    IF (omp_last) THEN
+      DO j=jj_begin-halo,jj_end+halo
+         DO i=ii_begin-halo,ii_end+halo
           dzq(ij,llm)=0.
-       ENDDO
-    ENDDO
+         ENDDO
+      ENDDO
+    ENDIF
+
+!---> flush dzq
+!$OMP BARRIER  
 
     ! sigw = fraction of mass that leaves level l/l+1
     ! then amount of q leaving level l/l+1 = wq = w * qq
-    DO l = 1,llm-1
+     DO l=ll_beginp1,ll_end
        DO j=jj_begin-halo,jj_end+halo
           DO i=ii_begin-halo,ii_end+halo
              ij=(j-1)*iim+i
-             w = fac*wfluxt(ij,l+1)
+             w = fac*wfluxt(ij,l)
              IF(w>0.) THEN  ! upward transport, upwind side is at level l 
+                sigw       = w/mass(ij,l-1)
+                qq         = q(ij,l-1)+0.5*(1.-sigw)*dzq(ij,l-1) ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0
+             ELSE           ! downward transport, upwind side is at level l+1
                 sigw       = w/mass(ij,l)
-                qq         = q(ij,l)+0.5*(1.-sigw)*dzq(ij,l) ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0
-             ELSE           ! downward transport, upwind side is at level l+1
-                sigw       = w/mass(ij,l+1)
-                qq         = q(ij,l+1)-0.5*(1.+sigw)*dzq(ij,l+1) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0                
+                qq         = q(ij,l)-0.5*(1.+sigw)*dzq(ij,l) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0                
              ENDIF
-             wq(ij,l+1) = w*qq
+             wq(ij,l) = w*qq
           ENDDO
        ENDDO
     END DO
     ! wq = 0 at top and bottom
-    DO j=jj_begin-halo,jj_end+halo
-       DO i=ii_begin-halo,ii_end+halo
-          ij=(j-1)*iim+i
-          wq(ij,llm+1)=0.
-          wq(ij,1)=0.
-       ENDDO
-    END DO
+    IF (omp_first) THEN
+      DO j=jj_begin-halo,jj_end+halo
+         DO i=ii_begin-halo,ii_end+halo
+            ij=(j-1)*iim+i
+            wq(ij,1)=0.
+         ENDDO
+      END DO
+    ENDIF
+    
+    IF (omp_last) THEN
+      DO j=jj_begin-halo,jj_end+halo
+         DO i=ii_begin-halo,ii_end+halo
+            ij=(j-1)*iim+i
+            wq(ij,llm+1)=0.
+         ENDDO
+      END DO
+    ENDIF
+
+! --> flush wq
+!$OMP BARRIER
+
 
     ! update q, mass is updated only after all q's have been updated
-    DO l=1,llm
+    DO l=ll_begin,ll_end
        DO j=jj_begin-halo,jj_end+halo
           DO i=ii_begin-halo,ii_end+halo

codes/icosagcm/trunk/src/caldyn.f90

@@ -11 +11 @@
   USE icosa
   USE caldyn_gcm_mod, ONLY : init_caldyn_gcm=>init_caldyn
-  USE caldyn_gcm_opt_mod, ONLY : init_caldyn_gcm_opt=>init_caldyn
   USE caldyn_adv_mod, ONLY : init_caldyn_adv=>init_caldyn
   IMPLICIT NONE
@@ -21 +20 @@
       CASE('gcm')
         CALL init_caldyn_gcm
-      CASE('gcm_opt')
-        CALL init_caldyn_gcm_opt
       CASE('adv')
         CALL init_caldyn_adv
@@ -37 +34 @@
   USE icosa
   USE caldyn_gcm_mod, ONLY : caldyn_gcm=>caldyn
-  USE caldyn_gcm_opt_mod, ONLY : caldyn_gcm_opt=>caldyn
   USE caldyn_adv_mod, ONLY : caldyn_adv=>caldyn
   IMPLICIT NONE
@@ -56 +52 @@
         CALL caldyn_gcm(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
              f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
-      CASE('gcm_opt')
-        CALL caldyn_gcm_opt(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
-             f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
       CASE('adv')
         CALL caldyn_adv(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &

codes/icosagcm/trunk/src/caldyn_gcm.f90

@@ -1 +1 @@
 MODULE caldyn_gcm_mod
   USE icosa
-
+  USE transfert_mod
   PRIVATE
 
@@ -11 +11 @@
   TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
 
-  PUBLIC init_caldyn, caldyn, write_output_fields,un2ulonlat
-
+! temporary shared variable for caldyn
+  TYPE(t_field),POINTER :: f_theta(:)
+  TYPE(t_field),POINTER :: f_pk(:)
+  TYPE(t_field),POINTER :: f_pks(:)
+  TYPE(t_field),POINTER :: f_phi(:)
+  TYPE(t_field),POINTER :: f_divm(:)
+  TYPE(t_field),POINTER :: f_wwuu(:)
+   
   INTEGER :: caldyn_hydrostat, caldyn_conserv
+  
+  TYPE(t_message) :: req_ps, req_theta_rhodz, req_u, req_qu
+
+  PUBLIC init_caldyn, caldyn, write_output_fields,req_ps
 
 CONTAINS
@@ -69 +79 @@
   IMPLICIT NONE
 
-  CALL allocate_field(f_out_u,field_u,type_real,llm) 
-  CALL allocate_field(f_p,field_t,type_real,llm+1)
-  CALL allocate_field(f_rhodz,field_t,type_real,llm) 
-  CALL allocate_field(f_qu,field_u,type_real,llm) 
-  
-  CALL allocate_field(f_buf_i,field_t,type_real,llm)
-  CALL allocate_field(f_buf_p,field_t,type_real,llm+1)
-  CALL allocate_field(f_buf_u3d,field_t,type_real,3,llm)  ! 3D vel at cell centers
-  CALL allocate_field(f_buf_ulon,field_t,type_real,llm)
-  CALL allocate_field(f_buf_ulat,field_t,type_real,llm)
-  CALL allocate_field(f_buf_v,field_z,type_real,llm)
-  CALL allocate_field(f_buf_s,field_t,type_real)
+    CALL allocate_field(f_out_u,field_u,type_real,llm) 
+    CALL allocate_field(f_p,field_t,type_real,llm+1)
+    CALL allocate_field(f_rhodz,field_t,type_real,llm) 
+    CALL allocate_field(f_qu,field_u,type_real,llm) 
+  
+    CALL allocate_field(f_buf_i,field_t,type_real,llm)
+    CALL allocate_field(f_buf_p,field_t,type_real,llm+1) 
+    CALL allocate_field(f_buf_u3d,field_t,type_real,3,llm)  ! 3D vel at cell centers
+    CALL allocate_field(f_buf_ulon,field_t,type_real,llm)
+    CALL allocate_field(f_buf_ulat,field_t,type_real,llm)
+    CALL allocate_field(f_buf_v,field_z,type_real,llm)
+    CALL allocate_field(f_buf_s,field_t,type_real)
+
+    CALL allocate_field(f_theta,field_t,type_real,llm)  ! potential temperature
+    CALL allocate_field(f_pk,field_t,type_real,llm)
+    CALL allocate_field(f_pks,field_t,type_real) ! Exner function
+    CALL allocate_field(f_phi,field_t,type_real,llm)    ! geopotential
+    CALL allocate_field(f_divm,field_t,type_real,llm)  ! mass flux divergence
+    CALL allocate_field(f_wwuu,field_u,type_real,llm+1) 
 
   END SUBROUTINE allocate_caldyn
@@ -92 +109 @@
     USE mpipara
     USE trace
+    USE omp_para
     IMPLICIT NONE
     LOGICAL,INTENT(IN)    :: write_out
@@ -108 +126 @@
     REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
     REAL(rstd),POINTER :: p(:,:), rhodz(:,:), qu(:,:)
+
+! temporary shared variable
+    REAL(rstd),POINTER  :: theta(:,:)  
+    REAL(rstd),POINTER  :: pk(:,:), pks(:)
+    REAL(rstd),POINTER  :: phi(:,:)    
+    REAL(rstd),POINTER  :: divm(:,:) 
+    REAL(rstd),POINTER  :: wwuu(:,:)
+    
+    
     INTEGER :: ind,ij
+    LOGICAL,SAVE :: first=.TRUE.
+!$OMP THREADPRIVATE(first)
+
+    
+    IF (first) THEN
+      first=.FALSE.
+      CALL init_message(f_ps,req_i1,req_ps)
+      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
+      CALL init_message(f_u,req_e1_vect,req_u)
+      CALL init_message(f_qu,req_e1_scal,req_qu)
+      CALL send_message(f_ps,req_ps) 
+    ENDIF
     
     CALL trace_start("caldyn")
-    CALL transfert_request(f_phis,req_i1) 
-    CALL transfert_request(f_ps,req_i1) 
-    CALL transfert_request(f_theta_rhodz,req_i1) 
-    CALL transfert_request(f_u,req_e1_vect)
+
+    CALL send_message(f_u,req_u)
+    CALL send_message(f_theta_rhodz,req_theta_rhodz) 
 
     SELECT CASE(caldyn_conserv)
@@ -126 +164 @@
           qu=f_qu(ind)
           u=f_u(ind)
-          !$OMP PARALLEL DEFAULT(SHARED)
+
           CALL compute_pvort(ps, u, p,rhodz,qu)
-          !$OMP END PARALLEL
+
        ENDDO
 
-       CALL transfert_request(f_qu,req_e1_scal)
+       CALL send_message(f_qu,req_qu)
 
        DO ind=1,ndomain
@@ -148 +186 @@
           u=f_u(ind)
           du=f_du(ind)
-          out_u=f_out_u(ind) 
-          !$OMP PARALLEL DEFAULT(SHARED)
-          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, &
-               hflux, wflux, dps, dtheta_rhodz, du)
-          !$OMP END PARALLEL
+          out_u=f_out_u(ind)
+          theta = f_theta(ind)
+          pk = f_pk(ind)
+          pks = f_pks(ind)
+          phi = f_phi(ind)  
+          divm = f_divm(ind)
+          wwuu=f_wwuu(ind)
+
+          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du, &
+                              theta,pk, pks, phi, divm, wwuu)
        ENDDO       
        
@@ -172 +215 @@
           du=f_du(ind)
           out_u=f_out_u(ind) 
-          !$OMP PARALLEL DEFAULT(SHARED)
+          wwuu=f_wwuu(ind)
           CALL compute_pvort(ps, u, p,rhodz,qu)
-          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, &
-               hflux, wflux, dps, dtheta_rhodz, du)
-          !$OMP END PARALLEL
+          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du,   &
+                              theta,pk, pks, phi, divm, wwuu)
        ENDDO
        
@@ -183 +225 @@
     END SELECT
 
+!$OMP BARRIER
     IF (write_out) THEN
+
        IF (is_mpi_root) PRINT *,'CALL write_output_fields'
-       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
-            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
+
+! ---> for openMP test to fix later
+!       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
+!            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
+       CALL writefield("ps",f_ps)
+       CALL writefield("dps",f_dps)
+       CALL writefield("vort",f_qu)
+       CALL writefield("theta",f_theta)
+
     END IF
     
     !    CALL check_mass_conservation(f_ps,f_dps)
     CALL trace_end("caldyn")
+!$OMP BARRIER
     
 END SUBROUTINE caldyn
@@ -199 +251 @@
   USE exner_mod
   USE trace
+  USE omp_para
   IMPLICIT NONE
   REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
@@ -208 +261 @@
   INTEGER :: i,j,ij,l
   REAL(rstd) :: etav,hv
-  REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)     ! potential velocity  
-  
-  LOGICAL,SAVE :: first=.TRUE.
-  !$OMP THREADPRIVATE(first)
-
+  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity  
+  
 
   CALL trace_start("compute_pvort")  
-  !$OMP BARRIER      
-  !$OMP MASTER  
-  !    IF (first) THEN
-  ALLOCATE(qv(2*iim*jjm,llm))     ! potential velocity  
-  
+
+  CALL wait_message(req_ps)   
 !!! Compute pressure
-  DO    l    = 1, llm+1
-     !$OMP DO
+   DO    l    = ll_begin, ll_endp1
+    CALL test_message(req_u) 
+
      DO j=jj_begin-1,jj_end+1
         DO i=ii_begin-1,ii_end+1
@@ -230 +278 @@
      ENDDO
   ENDDO
-  
+
+!$OMP BARRIER  
+
 !!! Compute mass
-  DO l = 1, llm
-     !$OMP DO
+   DO l = ll_begin,ll_end
+     CALL test_message(req_u) 
      DO j=jj_begin-1,jj_end+1
         DO i=ii_begin-1,ii_end+1
@@ -241 +291 @@
      ENDDO
   ENDDO
+
+  CALL wait_message(req_u)   
   
 !!! Compute shallow-water potential vorticity
-  DO l = 1,llm
-!$OMP DO
+  DO l = ll_begin,ll_end
+    CALL test_message(req_theta_rhodz) 
+
     DO j=jj_begin-1,jj_end+1
        DO i=ii_begin-1,ii_end+1
           ij=(j-1)*iim+i
           
-          etav= 1./Av(ij+z_up)*(  ne(ij,rup)        * u(ij+u_rup,l)        * de(ij+u_rup)         &
-                                + ne(ij+t_rup,left) * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
-                                - ne(ij,lup)        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
+          etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         &
+                                + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
+                                - ne_lup        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
 
           hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
@@ -259 +312 @@
           qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
           
-          etav = 1./Av(ij+z_down)*(  ne(ij,ldown)         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
-                                   + ne(ij+t_ldown,right) * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
-                                   - ne(ij,rdown)         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
+          etav = 1./Av(ij+z_down)*(  ne_ldown         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
+                                   + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
+                                   - ne_rdown         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
        
           hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
@@ -283 +336 @@
    ENDDO
 
-!!$OMP BARRIER
-!!$OMP MASTER  
-    DEALLOCATE(qv)       ! potential velocity  
-!!$OMP END MASTER
-!!$OMP BARRIER   
-    CALL trace_end("compute_pvort")
+   CALL trace_end("compute_pvort")
                                                        
   END SUBROUTINE compute_pvort
-  
-  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du)
+
+  
+  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du,   &
+                            theta,pk, pks, phi, divm, wwuu) 
   USE icosa
   USE disvert_mod
   USE exner_mod
   USE trace
+  USE omp_para
   IMPLICIT NONE
     REAL(rstd),INTENT(IN)  :: phis(iim*jjm)
@@ -310 +361 @@
     REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
     REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
-    
-    REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature
-    REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:), pks(:) ! Exner function
-    REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:), beta(:,:)
-    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential
-    REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux
-    REAL(rstd),ALLOCATABLE,SAVE :: divm(:,:)  ! mass flux divergence
-    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function
+
+! temporary variable
+    REAL(rstd),INTENT(INOUT) :: theta(iim*jjm,llm)  ! potential temperature
+    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm), pks(iim*jjm) ! Exner function
+    REAL(rstd),INTENT(INOUT) :: phi(iim*jjm,llm)    ! geopotential
+    REAL(rstd),INTENT(INOUT) :: divm(iim*jjm,llm)  ! mass flux divergence
+    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
+
+    
+    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
+    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernouilli function
     
     INTEGER :: i,j,ij,l
     REAL(rstd) :: ww,uu
 
-    LOGICAL,SAVE :: first=.TRUE.
-!$OMP THREADPRIVATE(first)
-    
     CALL trace_start("compute_caldyn")
-!$OMP BARRIER      
-!$OMP MASTER  
-!    IF (first) THEN
-    ALLOCATE(theta(iim*jjm,llm))    ! potential temperature
-    ALLOCATE(pk(iim*jjm,llm))       ! Exner function
-    ALLOCATE(pks(iim*jjm))
-    ALLOCATE(alpha(iim*jjm,llm))
-    ALLOCATE(beta(iim*jjm,llm))
-    ALLOCATE(phi(iim*jjm,llm))      ! geopotential
-    ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux   
-    ALLOCATE(divm(iim*jjm,llm))    ! mass flux divvergence
-    ALLOCATE(berni(iim*jjm,llm))    ! bernouilli term
+
+    CALL wait_message(req_theta_rhodz) 
 
 !!! Compute theta
-   DO l = 1, llm
-!$OMP DO
+   DO l = ll_begin, ll_end
+      IF (caldyn_conserv==energy) CALL test_message(req_qu) 
       DO j=jj_begin-1,jj_end+1
          DO i=ii_begin-1,ii_end+1
@@ -350 +391 @@
    ENDDO
 
-  DO l = 1, llm
+  DO l = ll_begin, ll_end
 !!!  Compute mass and theta fluxes
+    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
     DO j=jj_begin-1,jj_end+1
       DO i=ii_begin-1,ii_end+1
@@ -370 +412 @@
         ij=(j-1)*iim+i  
         ! divm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
-        divm(ij,l)= 1./Ai(ij)*(ne(ij,right)*hflux(ij+u_right,l)   +  &
-                                ne(ij,rup)*hflux(ij+u_rup,l)       +  &  
-                                ne(ij,lup)*hflux(ij+u_lup,l)       +  &  
-                                ne(ij,left)*hflux(ij+u_left,l)     +  &  
-                                ne(ij,ldown)*hflux(ij+u_ldown,l)   +  &  
-                                ne(ij,rdown)*hflux(ij+u_rdown,l))    
+        divm(ij,l)= 1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
+                                ne_rup*hflux(ij+u_rup,l)       +  &  
+                                ne_lup*hflux(ij+u_lup,l)       +  &  
+                                ne_left*hflux(ij+u_left,l)     +  &  
+                                ne_ldown*hflux(ij+u_ldown,l)   +  &  
+                                ne_rdown*hflux(ij+u_rdown,l))    
 
         ! signe ? attention d (rho theta dz)
         ! dtheta_rhodz =  -div(flux.theta)
-        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  &
-                                 ne(ij,rup)*Ftheta(ij+u_rup,l)        +  &  
-                                 ne(ij,lup)*Ftheta(ij+u_lup,l)        +  &  
-                                 ne(ij,left)*Ftheta(ij+u_left,l)      +  &  
-                                 ne(ij,ldown)*Ftheta(ij+u_ldown,l)    +  &  
-                                 ne(ij,rdown)*Ftheta(ij+u_rdown,l))
+        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  &
+                                 ne_rup*Ftheta(ij+u_rup,l)        +  &  
+                                 ne_lup*Ftheta(ij+u_lup,l)        +  &  
+                                 ne_left*Ftheta(ij+u_left,l)      +  &  
+                                 ne_ldown*Ftheta(ij+u_ldown,l)    +  &  
+                                 ne_rdown*Ftheta(ij+u_rdown,l))
       ENDDO
     ENDDO
   ENDDO
-   
+
+!$OMP BARRIER   
 !!! cumulate mass flux divergence from top to bottom
   DO  l = llm-1, 1, -1
-!$OMP DO
+    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
+!$OMP DO SCHEDULE(STATIC) 
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
@@ -399 +443 @@
     ENDDO
   ENDDO
+
+! IMPLICIT FLUSH on divm
+!!!!!!!!!!!!!!!!!!!!!!!!!
   
 !!! Compute vertical mass flux
-  DO l = 1,llm-1
-!$OMP DO
+!  DO l = 2,llm
+  DO l=ll_beginp1,ll_end
+    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
@@ -408 +456 @@
         ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
         ! => w>0 for upward transport
-        wflux( ij, l+1 ) = divm( ij, l+1 ) - bp(l+1) * divm( ij, 1 )
+        wflux( ij, l ) = divm( ij, l ) - bp(l) * divm( ij, 1 )
       ENDDO
     ENDDO
@@ -414 +462 @@
 
 ! compute dps, set vertical mass flux at the surface to 0
-!$OMP DO
-  DO j=jj_begin,jj_end
-    DO i=ii_begin,ii_end
-      ij=(j-1)*iim+i
-      wflux(ij,1)  = 0.
-      wflux(ij,llm+1)  = 0.
-      ! dps/dt = -int(div flux)dz
-      dps(ij)=-divm(ij,1) * g 
-    ENDDO
-  ENDDO
-
+  IF (omp_first) THEN
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        wflux(ij,1)  = 0.
+        ! dps/dt = -int(div flux)dz
+        dps(ij)=-divm(ij,1) * g 
+      ENDDO
+    ENDDO
+  ENDIF
+
+
+  IF (omp_last) THEN
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        wflux(ij,llm+1)  = 0.
+      ENDDO
+    ENDDO
+  ENDIF
 !!! Compute potential vorticity (Coriolis) contribution to du
 
@@ -430 +487 @@
     CASE(energy) ! energy-conserving TRiSK
 
-       DO l=1,llm
-          !$OMP DO
+       CALL wait_message(req_qu)
+
+        DO l=ll_begin,ll_end
           DO j=jj_begin,jj_end
              DO i=ii_begin,ii_end
@@ -479 +537 @@
     CASE(enstrophy) ! enstrophy-conserving TRiSK
   
-       DO l=1,llm
-          !$OMP DO
+        DO l=ll_begin,ll_end
           DO j=jj_begin,jj_end
              DO i=ii_begin,ii_end
@@ -531 +588 @@
   
 !!! Compute Exner function
-!    PRINT *, 'Computing Exner'
-    CALL compute_exner(ps,p,pks,pk,1) 
+!    CALL compute_exner(ps,p,pks,pk,1) 
+! replaced in source
+    IF (omp_first) THEN    
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+             ij=(j-1)*iim+i
+             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+          ENDDO
+       ENDDO
+     ENDIF
+     
+       ! 3D : pk
+     DO l = ll_begin, ll_end
+        DO j=jj_begin-1,jj_end+1
+           DO i=ii_begin-1,ii_end+1
+              ij=(j-1)*iim+i
+              pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
+           ENDDO
+        ENDDO
+     ENDDO
+
+!---> flush pk,pks, theta
+!$OMP BARRIER
 
 !!! Compute geopotential
 
   ! for first layer
-!$OMP DO
-   DO j=jj_begin-1,jj_end+1
-     DO i=ii_begin-1,ii_end+1
-       ij=(j-1)*iim+i
-       phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
+  IF (omp_first) THEN
+    DO j=jj_begin-1,jj_end+1
+      DO i=ii_begin-1,ii_end+1
+        ij=(j-1)*iim+i
+        phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
+      ENDDO
+    ENDDO
+  ENDIF
+!!-> implicit flush on phi(:,1)
+          
+  ! for other layers
+   DO l = ll_beginp1, ll_end
+     DO j=jj_begin-1,jj_end+1
+       DO i=ii_begin-1,ii_end+1
+         ij=(j-1)*iim+i
+         phi(ij,l) =  0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
+                          * (  pk(ij,l-1) -  pk(ij,l)    )
+       ENDDO
      ENDDO
    ENDDO
-          
-  ! for other layers
+
+!$OMP BARRIER
    DO l = 2, llm
 !$OMP DO
@@ -551 +642 @@
        DO i=ii_begin-1,ii_end+1
          ij=(j-1)*iim+i
-         phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
-                                       * (  pk(ij,l-1) -  pk(ij,l)    )
+         phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
        ENDDO
      ENDDO
    ENDDO
-
-   
+! --> IMPLICIT FLUSH on phi
+      
 !!! Compute bernouilli term = Kinetic Energy + geopotential
-  DO l=1,llm
-!$OMP DO
+   DO l=ll_begin,ll_end
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
@@ -576 +665 @@
     ENDDO
   ENDDO
-  
+
  
 !!! gradients of Bernoulli and Exner functions
-  DO l=1,llm
-!$OMP DO
+
+   DO l=ll_begin,ll_end
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
         ij=(j-1)*iim+i
         
-        out_u(ij+u_right,l)=  1/de(ij+u_right) * (                              &
-                          0.5*(theta(ij,l)+theta(ij+t_right,l))                             &
-                             *( ne(ij,right)*pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) &
-                        + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) )
+        du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             &
+                          0.5*(theta(ij,l)+theta(ij+t_right,l))                &
+                             *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    &
+                        + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
         
-        du(ij+u_right,l) = du(ij+u_right,l) + out_u(ij+u_right,l)
         
-        out_u(ij+u_lup,l)=  1/de(ij+u_lup) * (                                  &
-                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                             &
-                             *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l))    &
-                        + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) )
-        
-        du(ij+u_lup,l) = du(ij+u_lup,l) + out_u(ij+u_lup,l)
+        du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  &
+                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                  &
+                             *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       &
+                        + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
                 
-        out_u(ij+u_ldown,l)= 1/de(ij+u_ldown) * (                                &
-                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               &
-                             *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l))    &
-                        + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) )
-        
-        du(ij+u_ldown,l) = du(ij+u_ldown,l) + out_u(ij+u_ldown,l)
+        du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             &
+                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                &
+                             *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     &
+                        + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
 
       ENDDO
     ENDDO
   ENDDO
- 
-!!! contributions of vertical advection to du, dtheta
- 
-  DO l=1,llm-1
-!$OMP DO
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-         ! ww>0 <=> upward transport
-
-        ww            =   0.5 * wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
-        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww
-        dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww
-
-        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_right,l+1))
-        uu  = u(ij+u_right,l+1) - u(ij+u_right,l)  
-        du(ij+u_right, l )   = du(ij+u_right,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_right,l)))
-        du(ij+u_right, l+1 ) = du(ij+u_right,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_right,l+1)))
-
-        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_lup,l+1))
-        uu  = u(ij+u_lup,l+1) - u(ij+u_lup,l)  
-        du(ij+u_lup, l )   = du(ij+u_lup,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l)))
-        du(ij+u_lup, l+1 ) = du(ij+u_lup,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_lup,l+1)))
-
-        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_ldown,l+1))
-        uu  = u(ij+u_ldown,l+1) - u(ij+u_ldown,l)  
-        du(ij+u_ldown, l )   = du(ij+u_ldown,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l)))
-        du(ij+u_ldown, l+1 ) = du(ij+u_ldown,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_ldown,l+1)))
-
+
+   
+  DO l=ll_begin,ll_endm1
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
+      ENDDO
+    ENDDO
+  ENDDO
+  
+  DO l=ll_beginp1,ll_end
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) )
+      ENDDO
+    ENDDO
+  ENDDO
+  
+  IF (omp_first) THEN
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        wwuu(ij+u_right,1)=0   
+        wwuu(ij+u_lup,1)=0   
+        wwuu(ij+u_ldown,1)=0
+      ENDDO
+    ENDDO
+  ENDIF   
+
+  IF (omp_last) THEN
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        wwuu(ij+u_right,llm+1)=0   
+        wwuu(ij+u_lup,llm+1)=0   
+        wwuu(ij+u_ldown,llm+1)=0
+      ENDDO
+    ENDDO
+  ENDIF   
+  
+  DO l=ll_beginp1,ll_end
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        wwuu(ij+u_right,l) = 0.5*( wflux(ij,l) + wflux(ij+t_right,l)) * (u(ij+u_right,l) - u(ij+u_right,l-1))
+        wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1))
+        wwuu(ij+u_ldown,l) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1))
+       ENDDO
+     ENDDO
+   ENDDO
+
+ !--> flush wwuu  
+ !$OMP BARRIER
+
+  DO l=ll_begin,ll_end
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        du(ij+u_right, l )   = du(ij+u_right,l)  - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
+        du(ij+u_lup, l )     = du(ij+u_lup,l)    - (wwuu(ij+u_lup,l+1)  + wwuu(ij+u_lup,l))   / (rhodz(ij,l)+rhodz(ij+t_lup,l))
+        du(ij+u_ldown, l )   = du(ij+u_ldown,l)  - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
       ENDDO
     ENDDO
   ENDDO      
-   
-!!$OMP BARRIER
-!!$OMP MASTER  
-    DEALLOCATE(theta)   ! potential temperature
-    DEALLOCATE(pk)      ! Exner function
-    DEALLOCATE(pks)
-    DEALLOCATE(alpha)
-    DEALLOCATE(beta)
-    DEALLOCATE(phi)     ! geopotential
-    DEALLOCATE(Ftheta)  ! theta flux   
-    DEALLOCATE(divm)    ! mass flux divergence
-    DEALLOCATE(berni)   ! bernouilli term
-!!$OMP END MASTER
-!!$OMP BARRIER                                                      
-    CALL trace_end("compute_caldyn")
+  
+  CALL trace_end("compute_caldyn")
 
   END SUBROUTINE compute_caldyn
@@ -709 +814 @@
     USE write_field
     USE vertical_interp_mod
+    USE wind_mod
     TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
          f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
@@ -722 +828 @@
     
     CALL writefield("ps",f_ps)
-!    CALL writefield("dps",f_dps)
-!    CALL writefield("phis",f_phis)
-!    CALL vorticity(f_u,f_buf_v)
-!    CALL writefield("vort",f_buf_v)
-
-!    CALL w_omega(f_ps, f_u, f_buf_i)
-!    CALL writefield('omega', f_buf_i)
-!    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-!      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
-!      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
-!    ENDIF
+    CALL writefield("dps",f_dps)
+    CALL writefield("phis",f_phis)
+    CALL vorticity(f_u,f_buf_v)
+    CALL writefield("vort",f_buf_v)
+
+    CALL w_omega(f_ps, f_u, f_buf_i)
+    CALL writefield('omega', f_buf_i)
+    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
+      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
+      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
+    ENDIF
     
     ! Temperature
@@ -754 +860 @@
    
     ! velocity components
-    CALL un2ulonlat(f_u, f_buf_i3, f_buf1_i, f_buf2_i)
+    CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i)
     CALL writefield("ulon",f_buf1_i)
     CALL writefield("ulat",f_buf2_i)
@@ -767 +873 @@
     ! geopotential
     CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
-!    CALL writefield("p",f_buf_p)
-!    CALL writefield("phi",f_buf_i)
-     CALL writefield("theta",f_buf1_i) ! potential temperature
-!    CALL writefield("pk",f_buf2_i) ! Exner pressure
+    CALL writefield("p",f_buf_p)
+    CALL writefield("phi",f_buf_i)
+    CALL writefield("theta",f_buf1_i) ! potential temperature
+    CALL writefield("pk",f_buf2_i) ! Exner pressure
   
   
@@ -806 +912 @@
   END SUBROUTINE thetarhodz2geopot
   
-  SUBROUTINE un2ulonlat(f_u, f_u3d, f_ulon, f_ulat)
-    USE field_mod
-    USE wind_mod
-    TYPE(t_field), POINTER :: f_u(:), & ! IN  : normal velocity components on edges
-         f_u3d(:), f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons
-    REAL(rstd),POINTER :: u(:,:), u3d(:,:,:), ulon(:,:), ulat(:,:)
-    INTEGER :: ind
-    DO ind=1,ndomain
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       u=f_u(ind)
-       u3d=f_u3d(ind)
-       CALL compute_wind_centered(u,u3d)
-       ulon=f_ulon(ind)
-       ulat=f_ulat(ind)
-       CALL compute_wind_centered_lonlat_compound(u3d, ulon, ulat)
-    END DO
-  END SUBROUTINE un2ulonlat
-
   SUBROUTINE Tv2T(f_Tv, f_q, f_T)
   USE icosa

codes/icosagcm/trunk/src/check_conserve.f90

@@ -1 +1 @@
 MODULE check_conserve_mod
-        USE icosa 
-        IMPLICIT NONE 
-
-        PRIVATE 
-     TYPE(t_field),POINTER,SAVE :: f_pk(:), f_pks(:), f_p(:)
-     TYPE(t_field),POINTER,SAVE :: f_vort(:)
-     TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
+  USE icosa 
+  IMPLICIT NONE 
+
+  PRIVATE 
+    TYPE(t_field),POINTER,SAVE :: f_pk(:), f_pks(:), f_p(:)
+    TYPE(t_field),POINTER,SAVE :: f_vort(:)
+    TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
    
-     PUBLIC compute_conserve 
-     REAL(rstd),SAVE:: mtot0,ztot0,etot0,ang0,stot0,rmsv0
-     REAL(rstd),SAVE::rmsdpdt,etot,ang,stot,rmsv
-     REAL(rstd),SAVE:: ztot,mtot
+  PUBLIC init_check_conserve, check_conserve 
+    REAL(rstd),SAVE:: mtot0,ztot0,etot0,ang0,stot0,rmsv0
+    REAL(rstd),SAVE::rmsdpdt,etot,ang,stot,rmsv
+    REAL(rstd),SAVE:: ztot,mtot
         
      
-Contains 
-
-!---------------------------------------------------------------------
-
-   SUBROUTINE allocate_check_conserve
-    USE icosa
-    IMPLICIT NONE
+CONTAINS 
+
+!---------------------------------------------------------------------
+
+  SUBROUTINE init_check_conserve
+  USE icosa
+  IMPLICIT NONE
     CALL allocate_field(f_pk,field_t,type_real,llm)
     CALL allocate_field(f_p,field_t,type_real,llm+1)
@@ -26 +26 @@
     CALL allocate_field(f_rhodz,field_t,type_real,llm)
     CALL allocate_field(f_vort,field_z,type_real,llm)
-   END SUBROUTINE allocate_check_conserve
-
-!---------------------------------------------------------------------
-
-   SUBROUTINE check_mass_conserve(f_ps,f_dps)
+  END SUBROUTINE init_check_conserve
+
+  SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
+  USE icosa 
+  USE pression_mod
+  USE vorticity_mod
+  USE caldyn_gcm_mod
+  USE exner_mod 
+  USE mpipara
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_dps(:)
+    TYPE(t_field),POINTER :: f_ue(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_phis(:)
+    INTEGER::it 
+    REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
+    INTEGER::ind 
+
+    etot=0.0; ang=0.0;stot=0.0;rmsv=0.0 
+    mtot=0.0 ; rmsdpdt=0.0 ; ztot = 0.0 
+
+    CALL pression(f_ps,f_p)
+
+    DO ind=1,ndomain 
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+        p=f_p(ind) 
+       rhodz=f_rhodz(ind)
+      CALL compute_rhodz(p,rhodz) 
+    END DO 
+
+    CALL vorticity(f_ue,f_vort)
+    CALL check_mass_conserve(f_ps,f_dps)
+    CALL check_PV 
+    CALL exner(f_ps,f_p,f_pks,f_pk)
+    CALL check_EN(f_ue,f_theta_rhodz,f_phis) 
+       
+     IF ( it == 0  ) Then 
+       ztot0 = ztot
+       mtot0 = mtot
+       etot0 = etot 
+       ang0  = ang 
+       stot0 = stot 
+     END IF 
+
+     rmsv=SQRT(rmsv/mtot) 
+     ztot=ztot/ztot0 ; mtot=mtot/mtot0 
+     etot=etot/etot0 ; ang=ang/ang0 ; stot=stot/stot0 
+     rmsdpdt= daysec*1.e-2*sqrt(rmsdpdt/ncell_glo)  
+
+
+     IF (is_mpi_root) THEN 
+       OPEN(134,file="checkconsicosa.txt",position='append')
+       WRITE(134,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
+       WRITE(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang 
+       WRITE(134,*)"=================================================="
+       WRITE(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
+
+4000   FORMAT(10x,'masse',4x,'rmsdpdt',7x,'energie',2x,'enstrophie'  &
+     ,2x,'entropie',3x,'rmsv',4x,'mt.ang',/,'GLOB  '                &
+      ,f10.6,e13.6,5f10.3/)     
+       close(134) 
+     END IF 
+  END SUBROUTINE check_conserve
+  
+!---------------------------------------------------------------------
+
+  SUBROUTINE check_mass_conserve(f_ps,f_dps)
   USE icosa
   IMPLICIT NONE
@@ -38 +102 @@
     INTEGER :: ind,i,j,ij  
 
-   DO ind=1,ndomain
+    DO ind=1,ndomain
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
@@ -49 +113 @@
             mtot=mtot+ps(ij)*Ai(ij)
                rmsdpdt=rmsdpdt+dps(ij)*dps(ij)
-             ENDIF
+          ENDIF
         ENDDO
       ENDDO
@@ -57 +121 @@
 !---------------------------------------------------------------------
 
-   SUBROUTINE check_EN(f_ue,f_theta_rhodz,f_phis) 
-    USE icosa
-    USE pression_mod
-    USE vorticity_mod
-    IMPLICIT NONE
+  SUBROUTINE check_en(f_ue,f_theta_rhodz,f_phis) 
+  USE icosa
+  USE pression_mod
+  USE vorticity_mod
+  IMPLICIT NONE
     TYPE(t_field), POINTER :: f_ue(:)
     TYPE(t_field), POINTER :: f_theta_rhodz(:)
@@ -74 +138 @@
 
 
-      Do ind=1,ndomain 
-        CALL swap_dimensions(ind)
-        CALL swap_geometry(ind)
-         ue=f_ue(ind) 
-         p=f_p(ind)
-         rhodz=f_rhodz(ind) 
-         theta_rhodz=f_theta_rhodz(ind) 
-            pk=f_pk(ind) 
-            phis=f_phis(ind) 
-        CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis)
-         END DO 
-  END SUBROUTINE CHECK_EN 
-
-   SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis)
-         USE icosa
-      USE disvert_mod
-      USE wind_mod
-     IMPLICIT NONE
+    DO ind=1,ndomain 
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ue=f_ue(ind) 
+      p=f_p(ind)
+      rhodz=f_rhodz(ind) 
+      theta_rhodz=f_theta_rhodz(ind) 
+      pk=f_pk(ind) 
+      phis=f_phis(ind) 
+      CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis)
+    END DO 
+ 
+  END SUBROUTINE check_en 
+
+  SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis)
+  USE icosa
+  USE disvert_mod
+  USE wind_mod
+  IMPLICIT NONE
     INTEGER,INTENT(IN)::ind 
     REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm)
@@ -110 +175 @@
 
   
-       DO l = 1, llm
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-                ij=(j-1)*iim+i
-                 masse(ij,l) = rhodz(ij,l)*Ai(ij) 
-              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
-             IF (domain(ind)%own(i,j)) THEN
-                 stot = stot + Ai(ij)*theta_rhodz(ij,l)
-             END IF 
-          END DO 
-         END DO
-           END DO 
-
-       DO l = 1, llm
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-                ij=(j-1)*iim+i
-           IF (domain(ind)%own(i,j)) THEN
-                KE(ij,l)=   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   &
-                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           &
-                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
-                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
-                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
-                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
-             rmsv = rmsv + 2*masse(ij,l)*KE(ij,l) 
-            END IF 
-             END DO
-         END DO 
-        END DO  
-
-       DO l = 1, llm
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-                ij=(j-1)*iim+i
-           IF (domain(ind)%own(i,j)) THEN
-              etot = etot + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l))
-           END IF 
-          END DO 
+    DO l = 1, llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          masse(ij,l) = rhodz(ij,l)*Ai(ij) 
+          theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
+          IF (domain(ind)%own(i,j)) THEN
+            stot = stot + Ai(ij)*theta_rhodz(ij,l)
+          END IF 
         END DO 
-       END DO 
+      END DO
+    END DO 
+
+    DO l = 1, llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          IF (domain(ind)%own(i,j)) THEN
+            KE(ij,l)=   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   &
+                                      le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           &
+                                      le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
+                                      le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
+                                      le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
+                                      le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
+            rmsv = rmsv + 2*masse(ij,l)*KE(ij,l) 
+          END IF 
+        END DO
+      END DO 
+    END DO  
+
+    DO l = 1, llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          IF (domain(ind)%own(i,j)) THEN
+            etot = etot + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l))
+          END IF 
+        END DO 
+      END DO 
+    END DO 
 
 !------------------------------ ANGULAR VELOCITY 
-        CALL compute_wind_centered(u,ucenter)
-     CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
-         rad = radius 
-      radomeg = rad*omega
-
-      DO l = 1, llm
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-                ij=(j-1)*iim+i
-           IF (domain(ind)%own(i,j)) THEN
-             CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
-                ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat))
-           END IF 
-         END DO
-       END DO
-      END DO    
+    CALL compute_wind_centered(u,ucenter)
+    CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
+    rad = radius 
+    radomeg = rad*omega
+
+    DO l = 1, llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          IF (domain(ind)%own(i,j)) THEN
+            CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
+            ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat))
+          END IF 
+        END DO
+      END DO
+    END DO    
           
-        END SUBROUTINE compute_en
-
-!---------------------------------------------------------------------
-
-   SUBROUTINE check_PV 
-    USE icosa
-    IMPLICIT NONE
-
-    REAL(rstd), POINTER :: vort(:,:)
-    REAL(rstd), POINTER :: rhodz(:,:) 
-    INTEGER :: ind
-
-    Do ind=1,ndomain 
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      vort=f_vort(ind)
-      rhodz=f_rhodz(ind) 
-      CALL compute_PV(vort,rhodz) 
-    ENDDO
+  END SUBROUTINE compute_en
+
+!---------------------------------------------------------------------
+
+  SUBROUTINE check_PV 
+  USE icosa
+  IMPLICIT NONE
+     REAL(rstd), POINTER :: vort(:,:)
+     REAL(rstd), POINTER :: rhodz(:,:) 
+     INTEGER :: ind
+
+     DO ind=1,ndomain 
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       vort=f_vort(ind)
+       rhodz=f_rhodz(ind) 
+       CALL compute_PV(vort,rhodz) 
+     ENDDO
  
-   END SUBROUTINE check_PV 
-
-        SUBROUTINE compute_PV(vort,rhodz) 
-          USE icosa
-      USE disvert_mod
-     IMPLICIT NONE 
+  END SUBROUTINE check_PV 
+
+  SUBROUTINE compute_PV(vort,rhodz) 
+  USE icosa
+  USE disvert_mod
+  IMPLICIT NONE 
     REAL(rstd),INTENT(IN) :: vort(2*iim*jjm,llm)
     REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 
@@ -201 +265 @@
     INTEGER :: i,j,ij,l,ij2
 
-        hv1 = 0.0 ; hv2 = 0.0 
+    hv1 = 0.0 ; hv2 = 0.0 
    
-     DO l = 1,llm
+    DO l = 1,llm
       DO j=jj_begin+1,jj_end-1
-       DO i=ii_begin+1,ii_end-1
-           ij=(j-1)*iim+i
-            hv1 =  Riv2(ij,vup)        * rhodz(ij,l)           &
-              + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
-              + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
-      
-            qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
-            
-            hv2  =  Riv2(ij,vdown)     * rhodz(ij,l)          &
-              + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
-              + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
-                       
-            qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
-
-            ztot = ztot + (qv1*qv1*hv1 + qv2*qv2*hv2) 
-          ENDDO
-         ENDDO
+        DO i=ii_begin+1,ii_end-1
+          ij=(j-1)*iim+i
+          hv1 =  Riv2(ij,vup)        * rhodz(ij,l)           &
+            + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
+            + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
+    
+          qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
+          
+          hv2  =  Riv2(ij,vdown)     * rhodz(ij,l)          &
+            + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
+            + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
+                     
+          qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
+
+          ztot = ztot + (qv1*qv1*hv1 + qv2*qv2*hv2) 
+
         ENDDO
-    END SUBROUTINE compute_PV
-!---------------------------------------------------------------------
-
-    SUBROUTINE compute_rhodz(p,rhodz) 
-        USE icosa 
-        IMPLICIT NONE 
+      ENDDO
+    ENDDO
+
+  END SUBROUTINE compute_PV
+!---------------------------------------------------------------------
+
+   SUBROUTINE compute_rhodz(p,rhodz) 
+   USE icosa 
+   IMPLICIT NONE 
      REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
      REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm)
@@ -244 +310 @@
      ENDDO
 
-        END SUBROUTINE compute_rhodz
+   END SUBROUTINE compute_rhodz
 !====================================================================
 
-        SUBROUTINE compute_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
-         USE icosa 
-      USE pression_mod
-      USE vorticity_mod
-      USE caldyn_gcm_mod
-      USE exner_mod 
-      USE mpipara
-     TYPE(t_field),POINTER :: f_ps(:)
-     TYPE(t_field),POINTER :: f_dps(:)
-     TYPE(t_field),POINTER :: f_ue(:)
-     TYPE(t_field),POINTER :: f_theta_rhodz(:)
-     TYPE(t_field),POINTER :: f_phis(:)
-     INTEGER::it 
-     REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
-     INTEGER::ind 
-
-     etot=0.0; ang=0.0;stot=0.0;rmsv=0.0 
-     mtot=0.0 ; rmsdpdt=0.0 ; ztot = 0.0 
-
-     CALL allocate_check_conserve
-
-     CALL pression(f_ps,f_p)
-
-     Do ind=1,ndomain 
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-         p=f_p(ind) 
-        rhodz=f_rhodz(ind)
-       CALL compute_rhodz(p,rhodz) 
-      END DO 
-
-     CALL vorticity(f_ue,f_vort)
-     CALL check_mass_conserve(f_ps,f_dps)
-     CALL check_PV 
-     CALL exner(f_ps,f_p,f_pks,f_pk)
-     CALL check_EN(f_ue,f_theta_rhodz,f_phis) 
-        
-           IF ( it == 0  ) Then 
-           ztot0 = ztot
-           mtot0 = mtot
-           etot0 = etot 
-           ang0  = ang 
-           stot0 = stot 
-        END IF 
-
-      rmsv=SQRT(rmsv/mtot) 
-      ztot=ztot/ztot0 ; mtot=mtot/mtot0 
-      etot=etot/etot0 ; ang=ang/ang0 ; stot=stot/stot0 
-      rmsdpdt= daysec*1.e-2*sqrt(rmsdpdt/ncell_glo)  
-
-
-          IF (is_mpi_root) THEN 
-       OPEN(134,file="checkconsicosa.txt",position='append')
-       write(134,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
-       write(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang 
-       write(134,*)"=================================================="
-       write(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
-
-4000   FORMAT(10x,'masse',4x,'rmsdpdt',7x,'energie',2x,'enstrophie'  &
-      ,2x,'entropie',3x,'rmsv',4x,'mt.ang',/,'GLOB  '                &
-       ,f10.6,e13.6,5f10.3/)     
-       close(134) 
-       END IF 
-    END SUBROUTINE compute_conserve
 
  END MODULE check_conserve_mod

codes/icosagcm/trunk/src/dimensions.f90

@@ -44 +44 @@
     TYPE(t_domain),POINTER :: d
     
+    
+!$OMP MASTER
     d=>domain(ind)
      
@@ -79 +81 @@
      u_pos(:)=d%u_pos(:)
      z_pos(:)=d%z_pos(:)
-     
+
+!$OMP END MASTER     
+!$OMP BARRIER
    END SUBROUTINE swap_dimensions
    

codes/icosagcm/trunk/src/dissip_gcm.f90

@@ -10 +10 @@
   TYPE(t_field),POINTER,SAVE :: f_dtheta_diss(:)
   TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz_diss(:)
-  
+  TYPE(t_message) :: req_due, req_dtheta 
  
   INTEGER,SAVE :: nitergdiv=1
@@ -29 +29 @@
 !  INTEGER,SAVE :: itau_dissip
   REAL,SAVE    :: dtdissip
-  
+ 
   PUBLIC init_dissip, dissip
 CONTAINS
@@ -57 +57 @@
   USE mpi_mod
   USE mpipara
+  USE transfert_mod
   USE time_mod
   
@@ -96 +97 @@
       IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
    CASE DEFAULT
-      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), & 
-        ' in dissip_gcm.f90/init_dissip'
+      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), ' in dissip_gcm.f90/init_dissip'
       STOP
    END SELECT
@@ -116 +116 @@
     CALL allocate_field(f_theta,field_t,type_real)
     CALL allocate_field(f_dtheta,field_t,type_real)
+    
+    CALL init_message(f_due_diss1,req_e1_vect,req_due)
+    CALL init_message(f_dtheta_diss,req_i1,req_dtheta)
 
     tau_graddiv(:)=5000
@@ -135 +138 @@
 
     CALL getin("niterdivgrad",niterdivgrad)
-  
-!   CALL create_request(field_u,req_dissip)
-
-!   DO ind=1,ndomain
-!     DO i=ii_begin,ii_end
-!       CALL request_add_point(ind,i,jj_begin-1,req_dissip,rup)
-!       CALL request_add_point(ind,i+1,jj_begin-1,req_dissip,lup)
-!     ENDDO
-
-!     DO j=jj_begin,jj_end
-!       CALL request_add_point(ind,ii_end+1,j,req_dissip,left)
-!       CALL request_add_point(ind,ii_end+1,j-1,req_dissip,lup)
-!     ENDDO    
-!   
-!     DO i=ii_begin,ii_end
-!       CALL request_add_point(ind,i,jj_end+1,req_dissip,ldown)
-!       CALL request_add_point(ind,i-1,jj_end+1,req_dissip,rdown)
-!     ENDDO    
-
-!     DO j=jj_begin,jj_end
-!       CALL request_add_point(ind,ii_begin-1,j,req_dissip,right)
-!       CALL request_add_point(ind,ii_begin-1,j+1,req_dissip,rdown)
-!     ENDDO   
-
-!     DO i=ii_begin+1,ii_end-1
-!       CALL request_add_point(ind,i,jj_begin,req_dissip,right)
-!       CALL request_add_point(ind,i,jj_end,req_dissip,right)
-!     ENDDO
-!   
-!     DO j=jj_begin+1,jj_end-1
-!       CALL request_add_point(ind,ii_begin,j,req_dissip,rup)
-!       CALL request_add_point(ind,ii_end,j,req_dissip,rup)
-!     ENDDO   
-
-!     CALL request_add_point(ind,ii_begin+1,jj_begin,req_dissip,left)
-!     CALL request_add_point(ind,ii_begin,jj_begin+1,req_dissip,ldown)
-!     CALL request_add_point(ind,ii_begin+1,jj_end,req_dissip,left)
-!     CALL request_add_point(ind,ii_end,jj_begin+1,req_dissip,ldown)
-!    
-!   ENDDO
 
 
@@ -214 +177 @@
           u=f_u(ind)
           du=f_du(ind)
-          CALL compute_gradiv(u,du,1)
+          CALL compute_gradiv(u,du,1,1)
           u=du
         ENDDO
@@ -238 +201 @@
 
       IF (using_mpi) THEN 
-                          CALL MPI_ALLREDUCE(dumax,dumax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
+        CALL MPI_ALLREDUCE(dumax,dumax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
         dumax=dumax1
       ENDIF
@@ -291 +254 @@
           u=f_u(ind)
           du=f_du(ind)
-          CALL compute_gradrot(u,du,1)
+          CALL compute_gradrot(u,du,1,1)
           u=du
         ENDDO
@@ -362 +325 @@
           theta=f_theta(ind)
           dtheta=f_dtheta(ind)
-          CALL compute_divgrad(theta,dtheta,1)
+          CALL compute_divgrad(theta,dtheta,1,1)
           theta=dtheta
         ENDDO
       ENDDO
-!      CALL writefield("divgrad",f_dtheta)
 
       CALL transfert_request(f_dtheta,req_i1)
@@ -398 +360 @@
     ENDDO 
 
-!    CALL writefield("divgrad",f_dtheta)
     IF (is_mpi_root) PRINT *,"divgrad : divgrad",dthetamax
   
     cdivgrad=dthetamax**(-1./niterdivgrad) 
     IF (is_mpi_root) PRINT *,"cdivgrad : ",cdivgrad
-    
-
-
-
-
 
      
@@ -442 +398 @@
   USE trace
   USE time_mod
+  USE omp_para
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ue(:)
@@ -458 +415 @@
     INTEGER :: ind, shear
     INTEGER :: l,i,j,n
+
+!$OMP BARRIER
     
     CALL trace_start("dissip")
     CALL gradiv(f_ue,f_due_diss1)
     CALL gradrot(f_ue,f_due_diss2)
-    CALL theta_rhodz2theta(f_ps,f_theta_rhodz,f_theta)
-    CALL divgrad(f_theta,f_dtheta_diss)
-    CALL theta2theta_rhodz(f_ps,f_dtheta_diss,f_dtheta_rhodz_diss)
-    
-    IF(rayleigh_friction_type>0) THEN
-       CALL pression(f_ps, f_p)
-       CALL exner(f_ps, f_p, f_pks, f_pk)
-       CALL geopotential(f_phis,f_pks,f_pk,f_theta,f_phi)
-    END IF
+    
+    CALL divgrad_theta_rhodz(f_ps,f_theta_rhodz,f_dtheta_rhodz_diss)
+    
+! later for openmp    
+!    IF(rayleigh_friction_type>0) THEN
+!       CALL pression(f_ps, f_p)
+!       CALL exner(f_ps, f_p, f_pks, f_pk)
+!       CALL geopotential(f_phis,f_pks,f_pk,f_theta,f_phi)
+!    END IF
 
     DO ind=1,ndomain
@@ -480 +439 @@
       dtheta_rhodz=f_dtheta_rhodz(ind)
       dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind)
-            
-      DO l=1,llm
+
+      DO l=ll_begin,ll_end
         DO j=jj_begin,jj_end
           DO i=ii_begin,ii_end
@@ -495 +454 @@
       ENDDO
 
-      IF(rayleigh_friction_type>0) THEN
-         phi=f_phi(ind)
-         ue=f_ue(ind)
-         DO l=1,llm
-            DO j=jj_begin,jj_end
-               DO i=ii_begin,ii_end
-                  n=(j-1)*iim+i
-                  CALL relax(t_right, u_right)
-                  CALL relax(t_lup,   u_lup)
-                  CALL relax(t_ldown, u_ldown)
-               ENDDO
-            ENDDO
-         END DO
-      END IF
+!      dtheta_rhodz=0
+!      due=0
+
+! later for openmp  
+!      IF(rayleigh_friction_type>0) THEN
+!         phi=f_phi(ind)
+!         ue=f_ue(ind)
+!         DO l=1,llm
+!            DO j=jj_begin,jj_end
+!               DO i=ii_begin,ii_end
+!                  n=(j-1)*iim+i
+!                  CALL relax(t_right, u_right)
+!                  CALL relax(t_lup,   u_lup)
+!                  CALL relax(t_ldown, u_ldown)
+!               ENDDO
+!            ENDDO
+!         END DO
+!      END IF
    END DO
 
    CALL trace_end("dissip")
+
+!$OMP BARRIER
 
     CONTAINS
@@ -524 +489 @@
         z = (phi(n,l)+phi(n+shift_t,l))/(2.*g)
         IF(z>zh) THEN  ! relax only in the sponge layer z>zh
-!           PRINT *, 'z>zh : z,zh,l',z,zh,l
-!           STOP
+
            nn = n+shift_u
            CALL xyz2lonlat(xyz_e(nn,:),lon,lat)
@@ -548 +512 @@
   USE icosa
   USE trace
+  USE omp_para
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ue(:)
@@ -554 +519 @@
     REAL(rstd),POINTER    :: due(:,:)
     INTEGER :: ind
-    INTEGER :: it
+    INTEGER :: it,i,j,l,ij
        
     CALL trace_start("gradiv")
@@ -563 +528 @@
       ue=f_ue(ind)
       due=f_due(ind) 
-      due=ue
+      DO  l = ll_begin, ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+             due(ij+u_right,l)=ue(ij+u_right,l)
+             due(ij+u_lup,l)=ue(ij+u_lup,l)
+             due(ij+u_ldown,l)=ue(ij+u_ldown,l)
+          ENDDO
+        ENDDO
+      ENDDO
     ENDDO
 
     DO it=1,nitergdiv
         
-      CALL transfert_request(f_due,req_e1_vect)
+      CALL send_message(f_due,req_due)
+      CALL wait_message(req_due)
         
       DO ind=1,ndomain
@@ -574 +549 @@
         CALL swap_geometry(ind)
         due=f_due(ind) 
-        CALL compute_gradiv(due,due,llm)
+        CALL compute_gradiv(due,due,ll_begin,ll_end)
       ENDDO
     ENDDO
@@ -586 +561 @@
   USE icosa
   USE trace
+  USE omp_para
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ue(:)
@@ -592 +568 @@
     REAL(rstd),POINTER    :: due(:,:)
     INTEGER :: ind
-    INTEGER :: it
+    INTEGER :: it,i,j,l,ij
        
     CALL trace_start("gradrot")
@@ -601 +577 @@
       ue=f_ue(ind)
       due=f_due(ind) 
-      due=ue
+      DO  l = ll_begin, ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+             due(ij+u_right,l)=ue(ij+u_right,l)
+             due(ij+u_lup,l)=ue(ij+u_lup,l)
+             due(ij+u_ldown,l)=ue(ij+u_ldown,l)
+          ENDDO
+        ENDDO
+      ENDDO
     ENDDO
 
     DO it=1,nitergrot
         
-      CALL transfert_request(f_due,req_e1_vect)
+      CALL send_message(f_due,req_due)
+      CALL wait_message(req_due)
         
       DO ind=1,ndomain
@@ -612 +598 @@
         CALL swap_geometry(ind)
         due=f_due(ind) 
-        CALL compute_gradrot(due,due,llm)
+        CALL compute_gradrot(due,due,ll_begin,ll_end)
       ENDDO
 
@@ -624 +610 @@
   USE icosa
   USE trace
+  USE omp_para
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_theta(:)
@@ -650 +637 @@
         CALL swap_geometry(ind)
         dtheta=f_dtheta(ind) 
-        CALL compute_divgrad(dtheta,dtheta,llm)
+        CALL compute_divgrad(dtheta,dtheta,ll_begin,ll_end)
       ENDDO
 
@@ -659 +646 @@
   END SUBROUTINE divgrad
    
-  
-
-!  SUBROUTINE compute_dissip(ue,due,ps,theta_rhodz,dtheta_rhodz)
-!  USE icosa
-!  USE theta2theta_rhodz_mod
-!  IMPLICIT NONE
-!    REAL(rstd)         :: ue(3*iim*jjm,llm)
-!    REAL(rstd)         :: due(3*iim*jjm,llm)
-!    REAL(rstd)         :: ps(iim*jjm)
-!    REAL(rstd)         :: theta_rhodz(iim*jjm,llm)
-!    REAL(rstd)         :: dtheta_rhodz(iim*jjm,llm)
-!
-!    REAL(rstd),SAVE,ALLOCATABLE :: theta(:,:)
-!    REAL(rstd),SAVE,ALLOCATABLE :: du_dissip(:,:)
-!    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_dissip(:,:)
-!    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_rhodz_dissip(:,:)
-!
-!    INTEGER :: ind
-!    INTEGER :: l,i,j,n
-!
-!!$OMP BARRIER
-!!$OMP MASTER
-!      ALLOCATE(theta(iim*jjm,llm))
-!      ALLOCATE(du_dissip(3*iim*jjm,llm))
-!      ALLOCATE(dtheta_dissip(iim*jjm,llm))
-!      ALLOCATE(dtheta_rhodz_dissip(iim*jjm,llm))
-!!$OMP END MASTER
-!!$OMP BARRIER
-!    
-!      CALL gradiv(ue,du_dissip,llm)
-!      DO l=1,llm
-!!$OMP DO
-!        DO j=jj_begin,jj_end
-!          DO i=ii_begin,ii_end
-!            n=(j-1)*iim+i
-!            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_graddiv(l)*0.5
-!            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_graddiv(l)*0.5
-!            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_graddiv(l)*0.5
-!          ENDDO
-!        ENDDO
-!      ENDDO
-!      
-!      CALL gradrot(ue,du_dissip,llm)
-!
-!      DO l=1,llm
-!!$OMP DO
-!        DO j=jj_begin,jj_end
-!          DO i=ii_begin,ii_end
-!            n=(j-1)*iim+i
-!            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_gradrot(l)*0.5
-!            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_gradrot(l)*0.5
-!            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_gradrot(l)*0.5
-!          ENDDO
-!        ENDDO
-!      ENDDO 
-!      
-!      CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1)
-!      CALL divgrad(theta,dtheta_dissip,llm)
-!      CALL compute_theta2theta_rhodz(ps,dtheta_dissip,dtheta_rhodz_dissip,0)
-!
-!      DO l=1,llm
-!!$OMP DO
-!        DO j=jj_begin,jj_end
-!          DO i=ii_begin,ii_end
-!            n=(j-1)*iim+i
-!            dtheta_rhodz(n,l)=dtheta_rhodz(n,l)+dtheta_rhodz_dissip(n,l)/tau_divgrad(l)*0.5
-!          ENDDO
-!        ENDDO
-!      ENDDO 
-!
-!!$OMP BARRIER
-!!$OMP MASTER
-!      DEALLOCATE(theta)
-!      DEALLOCATE(du_dissip)
-!      DEALLOCATE(dtheta_dissip)
-!      DEALLOCATE(dtheta_rhodz_dissip)
-!!$OMP END MASTER
-!!$OMP BARRIER      
-!
-!  END SUBROUTINE compute_dissip
-  
-  
-  SUBROUTINE compute_gradiv(ue,gradivu_e,ll)
-  USE icosa
+  SUBROUTINE divgrad_theta_rhodz(f_ps,f_theta_rhodz,f_dtheta_rhodz)
+  USE icosa
+  USE trace
+  USE omp_para
+  USE disvert_mod
   IMPLICIT NONE
-    INTEGER,INTENT(IN)     :: ll
-    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,ll)
-    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,ll)
-    REAL(rstd) :: divu_i(iim*jjm,ll)
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
+    
+    REAL(rstd),POINTER :: ps(:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: dtheta_rhodz(:,:)
+
+    INTEGER :: ind
+    INTEGER :: it,i,j,l,ij
+
+    CALL trace_start("divgrad")
+       
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ps=f_ps(ind)
+      theta_rhodz=f_theta_rhodz(ind)
+      dtheta_rhodz=f_dtheta_rhodz(ind) 
+      DO  l = ll_begin, ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            dtheta_rhodz(ij,l) = theta_rhodz(ij,l) / ( ap(l)-ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g
+          ENDDO
+        ENDDO
+      ENDDO
+    ENDDO
+
+    DO it=1,niterdivgrad
+        
+      CALL send_message(f_dtheta_rhodz,req_dtheta)
+      CALL wait_message(req_dtheta)
+      DO ind=1,ndomain
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
+        dtheta_rhodz=f_dtheta_rhodz(ind) 
+        CALL compute_divgrad(dtheta_rhodz,dtheta_rhodz,ll_begin,ll_end)
+      ENDDO
+
+    ENDDO
+
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      dtheta_rhodz=f_dtheta_rhodz(ind) 
+    
+      DO  l = ll_begin, ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * ( ap(l)-ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g
+          ENDDO
+        ENDDO
+      ENDDO
+    ENDDO
+
+
+    CALL trace_end("divgrad")
+
+  END SUBROUTINE divgrad_theta_rhodz  
+  
+  SUBROUTINE compute_gradiv(ue,gradivu_e,llb,lle)
+  USE icosa
+  IMPLICIT NONE
+    INTEGER,INTENT(IN)     :: llb
+    INTEGER,INTENT(IN)     :: lle
+    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,llm)
+    REAL(rstd) :: divu_i(iim*jjm,llb:lle)
     
     INTEGER :: i,j,n,l
 
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -765 +739 @@
     ENDDO
     
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -781 +755 @@
     ENDDO
 
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -795 +769 @@
   END SUBROUTINE compute_gradiv
   
-  SUBROUTINE compute_divgrad(theta,divgrad_i,ll)
+  SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle)
   USE icosa
   IMPLICIT NONE
-    INTEGER,INTENT(IN)     :: ll
-    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,ll)
-    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,ll)
-    REAL(rstd)  :: grad_e(3*iim*jjm,ll)
+    INTEGER,INTENT(IN)     :: llb
+    INTEGER,INTENT(IN)     :: lle
+    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,llm)
+    REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle)
 
     INTEGER :: i,j,n,l
 
        
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin-1,jj_end+1
         DO i=ii_begin-1,ii_end+1
@@ -823 +798 @@
     
     
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -837 +812 @@
     ENDDO
    
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -849 +824 @@
 
     
-  SUBROUTINE compute_gradrot(ue,gradrot_e,ll)
+  SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle)
   USE icosa
   IMPLICIT NONE
-    INTEGER,INTENT(IN)     :: ll
-    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,ll)
-    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,ll)
-    REAL(rstd) :: rot_v(2*iim*jjm,ll)
+    INTEGER,INTENT(IN)     :: llb
+    INTEGER,INTENT(IN)     :: lle
+    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,llm)
+    REAL(rstd) :: rot_v(2*iim*jjm,llb:lle)
 
     INTEGER :: i,j,n,l
      
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin-1,jj_end+1
         DO i=ii_begin-1,ii_end+1
@@ -876 +852 @@
     ENDDO                              
   
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -891 +867 @@
     ENDDO
 
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end

codes/icosagcm/trunk/src/domain.f90

@@ -71 +71 @@
   SUBROUTINE create_domain
   USE grid_param
+  USE mpipara
   IMPLICIT NONE
   INTEGER :: ind,nf,ni,nj,i,j
@@ -98 +99 @@
           d%ii_end_glo=d%ii_begin_glo+d%ii_nb-1
  
-          IF (ni/=nsplit_i) THEN 
+          IF (ni/=1) THEN 
             d%ii_nb=d%ii_nb+1
-            d%ii_end_glo=d%ii_end_glo+1
+            d%ii_begin_glo=d%ii_begin_glo-1
           ENDIF
          
@@ -115 +116 @@
           d%jj_end_glo=d%jj_begin_glo+d%jj_nb-1
 
-          IF (nj/=nsplit_j) THEN 
+          IF (nj/=1) THEN 
             d%jj_nb=d%jj_nb+1
-            d%jj_end_glo=d%jj_end_glo+1
+            d%jj_begin_glo=d%jj_begin_glo-1
           ENDIF
 
@@ -142 +143 @@
           ALLOCATE(d%own(d%iim,d%jjm))
           ALLOCATE(d%ne(0:5,d%iim,d%jjm))
+          
+          IF (is_mpi_root) PRINT *,"Domain ",ind," : size ",d%ii_nb,"x",d%jj_nb
+          
         END DO
       END DO
@@ -261 +265 @@
             d%neighbour(:,k,i,j)=cell_glo(ind2)%xyz(:)
 
-!            d%ne(k,i,j)=vertex_glo(ii,jj,nf)%ne(k)
             d%ne(k,i,j)=1-2*MOD(k,2)
 
@@ -298 +301 @@
       edge_glo(e)%assign_pos=k
       edge_glo(e)%assign_delta=delta
-!      PRINT*,"Assign_edge",ind_d,ind,i,j,k,e
+
      END  SUBROUTINE assign_edge
          

codes/icosagcm/trunk/src/dynetat0_gcm_mod.f90

@@ -18 +18 @@
 CONTAINS
 
-        SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
+  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
   USE icosa
-  USE caldyn_mod
   USE write_field  
   USE maxicosa
@@ -98 +97 @@
 !----------------------------------------------------
 !------------- OUTPUT OF Variables 
-        CALL un2ulonlat(f_u,f_buf_i3,f_buf1_i,f_buf2_i) 
-        CALL writefield("buf1",f_buf1_i)
   END SUBROUTINE etat0
 

codes/icosagcm/trunk/src/dynetat0_hz_mod.f90

@@ -2 +2 @@
   USE genmod
   USE icosa
-  USE caldyn_gcm_mod 
   IMPLICIT NONE 
   PRIVATE
@@ -23 +22 @@
   USE write_field  
   USE maxicosa
+  USE wind_mod
         IMPLICIT NONE
      TYPE(t_domain),POINTER :: d 
@@ -97 +97 @@
         ENDDO
 !------------- OUTPUT OF Variables 
-        CALL un2ulonlat(f_u,f_buf_i3,f_buf1_i,f_buf2_i) 
+        CALL un2ulonlat(f_u,f_buf1_i,f_buf2_i) 
         CALL writefield("buf1",f_buf1_i)
   END SUBROUTINE etat0

codes/icosagcm/trunk/src/exner.f90

@@ -50 +50 @@
        
        ! for llm layer
-       !$OMP DO
        DO j=jj_begin-offset,jj_end+offset
           DO i=ii_begin-offset,ii_end+offset
@@ -61 +60 @@
        ! for other layer   
        DO l = llm-1 , 2 , -1
-          !$OMP DO
           DO j=jj_begin-offset,jj_end+offset
              DO i=ii_begin-offset,ii_end+offset
@@ -75 +73 @@
        
        ! for first layer
-       !$OMP DO
        DO j=jj_begin-offset,jj_end+offset
           DO i=ii_begin-offset,ii_end+offset
@@ -88 +85 @@
        
        DO l = 2, llm
-          !$OMP DO
           DO j=jj_begin-offset,jj_end+offset
              DO i=ii_begin-offset,ii_end+offset
@@ -99 +95 @@
     ELSE ! Simple calculation of Exner pressure based on centered average
        ! surface : pks
-       !$OMP DO
        DO j=jj_begin-offset,jj_end+offset
           DO i=ii_begin-offset,ii_end+offset
@@ -109 +104 @@
        ! 3D : pk
        DO l = 1, llm
-          !$OMP DO
           DO j=jj_begin-offset,jj_end+offset
              DO i=ii_begin-offset,ii_end+offset

codes/icosagcm/trunk/src/geometry.f90

@@ -91 +91 @@
   IMPLICIT NONE
     INTEGER,INTENT(IN) :: ind
+!$OMP MASTER
     Ai=geom%Ai(ind)
     xyz_i=geom%xyz_i(ind)
@@ -111 +112 @@
     bi=geom%bi(ind)
     fv=geom%fv(ind)
-    
+!$OMP END MASTER
+!$OMP BARRIER    
   END SUBROUTINE swap_geometry
   
@@ -213 +215 @@
           ENDDO
         ENDDO
-!        i=ii_begin ; j=jj_begin ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j) 
-!        i=ii_begin ; j=jj_end   ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j)
-!        i=ii_end   ; j=jj_begin ; n=(j-1)*iim+i ;   xyz_i(n,:)=domain(ind)%xyz(:,i,j)
-!        PRINT *,"Pb ?? : " 
-!        PRINT *, xyz_i(n,:), domain(ind)%xyz(:,i,j), norm(xyz_i(n,:)- domain(ind)%xyz(:,i,j))
-!        i=ii_end   ; j=jj_end   ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j)
         
       ENDDO
 
       IF (check) THEN
-!        sum=sum/(ndomain*ii_nb*jj_nb)
         PRINT *,"it = ",it,"  diff centroid circumcenter ",sum
       ENDIF 
@@ -287 +282 @@
             ne(n,k+1)=d%ne(k,i,j)
           ENDDO
-                  
-!          xyz_i(n,:)=d%xyz(:,i,j) 
-!          xyz_v(n+z_up,:)=d%vertex(:,vup-1,i,j) 
-!          xyz_v(n+z_down,:)=d%vertex(:,vdown-1,i,j) 
 
           vect(:,1)=xyz_v(n+z_rup,:)
@@ -309 +300 @@
           bi(n)=0. 
         
-!          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),de(n+u_right))
-!          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),de(n+u_lup))
-!          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),de(n+u_ldown))
-
           CALL dist_cart(xyz_i(n,:),xyz_i(n+t_right,:),de(n+u_right))
           CALL dist_cart(xyz_i(n,:),xyz_i(n+t_lup,:),de(n+u_lup))
           CALL dist_cart(xyz_i(n,:),xyz_i(n+t_ldown,:),de(n+u_ldown))
           
-!          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),0.5,xyz_e(n+u_right,:))
-!          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),0.5,xyz_e(n+u_lup,:))
-!          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),0.5,xyz_e(n+u_ldown,:))
-
           CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_right,:),0.5,xyz_e(n+u_right,:))
           CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_lup,:),0.5,xyz_e(n+u_lup,:))
           CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_ldown,:),0.5,xyz_e(n+u_ldown,:))
-         
-!          CALL dist_cart(d%vertex(:,vrdown-1,i,j),d%vertex(:,vrup-1,i,j),le(n+u_right))
-!          CALL dist_cart(d%vertex(:,vup-1,i,j),d%vertex(:,vlup-1,i,j),le(n+u_lup))
-!          CALL dist_cart(d%vertex(:,vldown-1,i,j),d%vertex(:,vdown-1,i,j),le(n+u_ldown))
 
           CALL dist_cart(xyz_v(n+z_rdown,:), xyz_v(n+z_rup,:),le(n+u_right))
@@ -335 +314 @@
           Ai(n)=0
           DO k=0,5
-!            CALL surf_triangle(d%xyz(:,i,j),d%neighbour(:,k,i,j),d%neighbour(:,MOD((k+1+6),6),i,j),surf_v(k+1))
-!            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,MOD((k-1+6),6),i,j),d%vertex(:,k,i,j),surf(k+1))
             CALL surf_triangle(xyz_i(n,:),xyz_i(n+t_pos(k+1),:),xyz_i(n+t_pos(MOD((k+1+6),6)+1),:),surf_v(k+1))
             CALL surf_triangle(xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:),surf(k+1))
@@ -414 +391 @@
           
           DO k=0,5
-!            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,k,i,j),S1)
-!            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,MOD(k+1+6,6),i,j),S2)
             CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1)
             CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(MOD(k+1+6,6)+1),:),S2)
@@ -516 +491 @@
     CALL transfert_request(geom%centroid,req_i1)
     CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S)
-!    PRINT *,"Surf triangle : ",S*20/(4*Pi)
-  
+ 
   END SUBROUTINE set_geometry
   

codes/icosagcm/trunk/src/icosa_gcm.f90

@@ -6 +6 @@
   USE wind_mod
   USE mpipara
+  USE omp_para
   USE vertical_interp_mod
+  USE trace
   IMPLICIT NONE
   
@@ -21 +23 @@
   CALL init_earth_const 
   CALL init_grid_param
+  CALL init_omp_para
   CALL compute_metric
   CALL compute_domain
   CALL init_transfert
   CALL init_writefield
+  CALL init_trace
 
-!  CALL allocate_field(sum_ne,field_T,type_real)
-!  CALL allocate_field_glo(sum_ne_glo,field_T,type_real)
-!  
-! DO ind=1,ndomain
-!   CALL swap_dimensions(ind)
-!   pt_sum_ne=sum_ne(ind)
-!   DO j=jj_begin,jj_end
-!     DO i=ii_begin,ii_end    
-!       n=(j-1)*iim+i
-!       pt_sum_ne(n)=domloc_glo_ind(ind)
-!     ENDDO
-!   ENDDO
-! ENDDO
-
-! CALL WriteField("domain",sum_ne)
-! CALL WriteField_mpi("domain",sum_ne)
-! CALL transfert_request(sum_ne,req_i1)
-! CALL WriteField_mpi("domain",sum_ne)
-! CALL close_files
-! CALL finalize_mpipara
-! STOP
   
   CALL compute_geometry
@@ -80 +63 @@
   
   CALL WriteField("Ai",geom%Ai)
-!  CALL WriteField("sum_ne",sum_ne)
+
   IF (is_mpi_root) CALL write_apbp
   CALL init_time
+
+  CALL init_timeloop
+
+!$OMP PARALLEL
   CALL timeloop
+!$OMP END PARALLEL
 
   CALL close_files

codes/icosagcm/trunk/src/mpi_mod.F90

@@ -11 +11 @@
   INTEGER :: MPI_INFO_NULL
   INTEGER :: MPI_STATUS_SIZE 
+  INTEGER,PARAMETER :: MPI_ADDRESS_KIND=KIND(INTEGER)
 #endif
 
@@ -40 +41 @@
  END
 
+ SUBROUTINE MPI_ISSEND
+ END
+
  SUBROUTINE MPI_IRECV
  END
 
  SUBROUTINE MPI_WAITALL
+ END
+
+ SUBROUTINE MPI_TESTALL
  END
 
@@ -51 +58 @@
  SUBROUTINE MPI_ALLGATHER
  END
+ 
+ SUBROUTINE MPI_TYPE_EXTENT
+ END
+ 
+ SUBROUTINE MPI_ALLOC_MEM
+ END
 
 #endif

codes/icosagcm/trunk/src/mpipara.F90

@@ -8 +8 @@
   LOGICAL,SAVE :: using_mpi
   LOGICAL,SAVE :: is_mpi_root
+  
+  INTERFACE allocate_mpi_buffer
+    MODULE PROCEDURE allocate_mpi_buffer_r2, allocate_mpi_buffer_r3,allocate_mpi_buffer_r4
+  END INTERFACE allocate_mpi_buffer
 
 CONTAINS
@@ -48 +52 @@
    END SUBROUTINE  finalize_mpipara
    
+
+  SUBROUTINE allocate_mpi_buffer_r2(buffer,length)
+  USE ISO_C_BINDING
+  USE mpi_mod
+  USE prec
+  IMPLICIT NONE
+    REAL(rstd), POINTER :: buffer(:)
+    INTEGER,INTENT(IN)  :: length
+
+    TYPE(C_PTR)         :: base_ptr
+    INTEGER(KIND=MPI_ADDRESS_KIND) :: size
+    INTEGER :: real_size,ierr
+    
+    CALL MPI_TYPE_EXTENT(MPI_REAL8, real_size, ierr)
+    size=length*real_size
+    
+    CALL MPI_ALLOC_MEM(size,MPI_INFO_NULL,base_ptr,ierr)
+    CALL C_F_POINTER(base_ptr, buffer, (/ length /))
+
+   END SUBROUTINE allocate_mpi_buffer_r2
+
+  SUBROUTINE allocate_mpi_buffer_r3(buffer,length,dim3)
+  USE ISO_C_BINDING
+  USE mpi_mod
+  USE prec
+    IMPLICIT NONE
+    REAL(rstd), POINTER :: buffer(:,:)
+    INTEGER,INTENT(IN)  :: length
+    INTEGER,INTENT(IN)  :: dim3
+
+    TYPE(C_PTR)         :: base_ptr
+    INTEGER(KIND=MPI_ADDRESS_KIND) :: size
+    INTEGER :: real_size,ierr
+    
+    CALL MPI_TYPE_EXTENT(MPI_REAL8, real_size, ierr)
+    size=length*real_size*dim3
+    
+    CALL MPI_ALLOC_MEM(size,MPI_INFO_NULL,base_ptr,ierr)
+    CALL C_F_POINTER(base_ptr, buffer, (/ length,dim3 /))
+    
+   END SUBROUTINE allocate_mpi_buffer_r3
+
+  SUBROUTINE allocate_mpi_buffer_r4(buffer,length,dim3,dim4)
+  USE ISO_C_BINDING
+  USE mpi_mod
+  USE prec
+  IMPLICIT NONE
+    REAL(rstd), POINTER :: buffer(:,:,:)
+    INTEGER,INTENT(IN)  :: length
+    INTEGER,INTENT(IN)  :: dim3
+    INTEGER,INTENT(IN)  :: dim4
+
+    TYPE(C_PTR)         :: base_ptr
+    INTEGER(KIND=MPI_ADDRESS_KIND) :: size
+    INTEGER :: real_size,ierr
+    
+    CALL MPI_TYPE_EXTENT(MPI_REAL8, real_size, ierr)
+    size=length*real_size*dim3*dim4
+    
+    CALL MPI_ALLOC_MEM(size,MPI_INFO_NULL,base_ptr,ierr)
+    CALL C_F_POINTER(base_ptr, buffer, (/ length, dim3, dim4 /))
+    
+   END SUBROUTINE allocate_mpi_buffer_r4
    
 END MODULE mpipara

codes/icosagcm/trunk/src/pression.f90

@@ -33 +33 @@
 
     DO    l    = 1, llm+1
-!$OMP DO
       DO j=jj_begin-offset,jj_end+offset
         DO i=ii_begin-offset,ii_end+offset

codes/icosagcm/trunk/src/theta_rhodz.f90

@@ -83 +83 @@
     REAL(rstd),ALLOCATABLE,SAVE :: p(:,:)
 
-!$OMP BARRIER
-!$OMP MASTER
-      ALLOCATE( p(iim*jjm,llm+1))
-!$OMP END MASTER
-!$OMP BARRIER
-
+    ALLOCATE( p(iim*jjm,llm+1))
     CALL compute_pression(ps,p,offset)
     
     DO    l    = 1, llm
-!$OMP DO
       DO j=jj_begin-offset,jj_end+offset
         DO i=ii_begin-offset,ii_end+offset
@@ -101 +95 @@
     ENDDO
 
-!$OMP BARRIER
-!$OMP MASTER
-      DEALLOCATE( p)
-!$OMP END MASTER
-!$OMP BARRIER    
+    DEALLOCATE( p)
 
   END SUBROUTINE compute_theta2theta_rhodz
@@ -120 +110 @@
     REAL(rstd),SAVE,ALLOCATABLE :: p(:,:)
 
-!$OMP BARRIER
-!$OMP MASTER
-      ALLOCATE( p(iim*jjm,llm+1))
-!$OMP END MASTER
-!$OMP BARRIER
+    ALLOCATE( p(iim*jjm,llm+1))
+
     CALL compute_pression(ps,p,offset)
     
@@ -136 +123 @@
     ENDDO
 
-!$OMP BARRIER
-!$OMP MASTER
-      DEALLOCATE( p)
-!$OMP END MASTER
-!$OMP BARRIER    
+    DEALLOCATE( p)
     
   END SUBROUTINE compute_theta_rhodz2theta
@@ -193 +176 @@
         DO i=ii_begin-offset,ii_end+offset
           ij=(j-1)*iim+i
-!          temp(ij,l) = ( theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g) ) * pk(ij,l) / cpp
           theta_rhodz(ij,l) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp )
         ENDDO

codes/icosagcm/trunk/src/time.f90

@@ -119 +119 @@
   INTEGER :: status
   REAL(rstd) ::time_array(1)
+
+!$OMP MASTER
     time_array(1)=time
   
@@ -126 +128 @@
       status=NF90_SYNC(ncid)
     ENDIF
+!$OMP END MASTER
 
   END SUBROUTINE update_time_counter    

codes/icosagcm/trunk/src/timeloop_gcm.f90

@@ -1 +1 @@
 MODULE timeloop_gcm_mod
-
+  USE transfert_mod
+  USE icosa
   PRIVATE
   
-  PUBLIC :: timeloop
+  PUBLIC :: init_timeloop, timeloop
 
   INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3
   INTEGER  :: itau_sync=10
 
+  TYPE(t_message) :: req_ps0, req_theta_rhodz0, req_u0, req_q0
+
+  TYPE(t_field),POINTER :: f_phis(:)
+  TYPE(t_field),POINTER :: f_q(:)
+  TYPE(t_field),POINTER :: f_dtheta(:), f_rhodz(:)
+  TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:)
+  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:)
+  TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:)
+  TYPE(t_field),POINTER :: f_dps(:),f_dpsm1(:), f_dpsm2(:)
+  TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)
+  TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:)
+  TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)  
+
+  INTEGER :: nb_stage, matsuno_period, scheme    
+
+  REAL(rstd),SAVE :: jD_cur, jH_cur
+  REAL(rstd),SAVE :: start_time
+
 CONTAINS
   
-  SUBROUTINE timeloop
+  SUBROUTINE init_timeloop
   USE icosa
   USE dissip_gcm_mod
@@ -19 +38 @@
   USE physics_mod
   USE mpipara
-  USE IOIPSL
-  USE maxicosa
-  USE check_conserve_mod
+  USE omp_para
   USE trace
   USE transfert_mod
+  USE check_conserve_mod
+  USE ioipsl
   IMPLICIT NONE
-  TYPE(t_field),POINTER :: f_phis(:)
-!  TYPE(t_field),POINTER :: f_theta(:)
-  TYPE(t_field),POINTER :: f_q(:)
-  TYPE(t_field),POINTER :: f_dtheta(:), f_rhodz(:)
-  TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:)
-  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:)
-  TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:)
-  TYPE(t_field),POINTER :: f_dps(:),f_dpsm1(:), f_dpsm2(:)
-  TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)
-  TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:)
-  TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)
-
-  REAL(rstd),POINTER :: phis(:)
-  REAL(rstd),POINTER :: q(:,:,:)
-  REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:)
-  REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:)
-  REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
-  REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:)
-  REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
-  REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
-  REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
-!  REAL(rstd) :: dt, run_length
-  INTEGER :: ind
-  INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme
-  CHARACTER(len=255) :: scheme_name
-  LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
-  CHARACTER(len=7) :: first 
-  REAL(rstd),SAVE :: jD_cur, jH_cur
-  REAL(rstd),SAVE :: start_time
-  LOGICAL, PARAMETER :: check=.FALSE.
-!  INTEGER :: itaumax
-!  REAL(rstd) ::write_period
-!  INTEGER    :: itau_out
-  
-!  dt=90.
-!  CALL getin('dt',dt)
-  
-!  itaumax=100
-!  CALL getin('itaumax',itaumax)
-
-!  run_length=dt*itaumax
-!  CALL getin('run_length',run_length)
-!  itaumax=run_length/dt
-!  PRINT *,'itaumax=',itaumax
-!  dt=dt/scale_factor
-
-! Trends
-  CALL allocate_field(f_dps,field_t,type_real)
-  CALL allocate_field(f_du,field_u,type_real,llm)
-  CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm)
+
+    CHARACTER(len=255) :: scheme_name
+
+    CALL allocate_field(f_dps,field_t,type_real)
+    CALL allocate_field(f_du,field_u,type_real,llm)
+    CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm)
 ! Model state at current time step (RK/MLF/Euler)
-  CALL allocate_field(f_phis,field_t,type_real)
-  CALL allocate_field(f_ps,field_t,type_real)
-  CALL allocate_field(f_u,field_u,type_real,llm)
-  CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
+    CALL allocate_field(f_phis,field_t,type_real)
+    CALL allocate_field(f_ps,field_t,type_real)
+    CALL allocate_field(f_u,field_u,type_real,llm)
+    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
 ! Model state at previous time step (RK/MLF)
-  CALL allocate_field(f_psm1,field_t,type_real)
-  CALL allocate_field(f_um1,field_u,type_real,llm)
-  CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm)
+    CALL allocate_field(f_psm1,field_t,type_real)
+    CALL allocate_field(f_um1,field_u,type_real,llm)
+    CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm)
 ! Tracers
-  CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
-  CALL allocate_field(f_rhodz,field_t,type_real,llm)
+    CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
+    CALL allocate_field(f_rhodz,field_t,type_real,llm)
 ! Mass fluxes
-  CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
-  CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! computed by caldyn
-  CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
-  CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
+    CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
+    CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! computed by caldyn
+    CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
+    CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
+
 
 !----------------------------------------------------
@@ -145 +121 @@
 !----------------------------------------------------
 
-  scheme_name='runge_kutta'
-  CALL getin('scheme',scheme_name)
-
-  SELECT CASE (TRIM(scheme_name))
-  CASE('euler')
-     scheme=euler
-     nb_stage=1
-  CASE ('runge_kutta')
-     scheme=rk4
-     nb_stage=4
-  CASE ('leapfrog_matsuno')
-     scheme=mlf
-     matsuno_period=5
-     CALL getin('matsuno_period',matsuno_period)
-     nb_stage=matsuno_period+1
+
+
+
+
+    scheme_name='runge_kutta'
+    CALL getin('scheme',scheme_name)
+
+    SELECT CASE (TRIM(scheme_name))
+      CASE('euler')
+        scheme=euler
+        nb_stage=1
+      CASE ('runge_kutta')
+        scheme=rk4
+        nb_stage=4
+      CASE ('leapfrog_matsuno')
+        scheme=mlf
+        matsuno_period=5
+        CALL getin('matsuno_period',matsuno_period)
+        nb_stage=matsuno_period+1
      ! Model state 2 time steps ago (MLF)
-     CALL allocate_field(f_psm2,field_t,type_real)
-     CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
-     CALL allocate_field(f_um2,field_u,type_real,llm)
-  CASE default
-     PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name),             &
-          ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
-     STOP
-  END SELECT
+        CALL allocate_field(f_psm2,field_t,type_real)
+        CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
+        CALL allocate_field(f_um2,field_u,type_real,llm)
+    CASE default
+       PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name),             &
+            ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
+       STOP
+    END SELECT
+
+
+    CALL init_dissip
+    CALL init_caldyn
+    CALL init_guided
+    CALL init_advect_tracer
+    CALL init_check_conserve
+    CALL init_physics
+    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+
+    CALL transfert_request(f_phis,req_i0) 
+    CALL transfert_request(f_phis,req_i1) 
+    CALL writefield("phis",f_phis,once=.TRUE.)
+
+    CALL init_message(f_ps,req_i0,req_ps0)
+    CALL init_message(f_theta_rhodz,req_i0,req_theta_rhodz0)
+    CALL init_message(f_u,req_e0_vect,req_u0)
+    CALL init_message(f_q,req_i0,req_q0)
+
+  END SUBROUTINE init_timeloop
   
-!  write_period=0
-!  CALL getin('write_period',write_period)
-!  write_period=write_period/scale_factor
-!  itau_out=FLOOR(.5+write_period/dt)
-!  PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out
- 
-! Trends at previous time steps needed only by Adams-Bashforth
-!  CALL allocate_field(f_dpsm1,field_t,type_real)
-!  CALL allocate_field(f_dpsm2,field_t,type_real)
-!  CALL allocate_field(f_dum1,field_u,type_real,llm)
-!  CALL allocate_field(f_dum2,field_u,type_real,llm)
-!  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm)
-!  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm)
-!  CALL allocate_field(f_theta,field_t,type_real,llm)
-!  CALL allocate_field(f_dtheta,field_t,type_real,llm)
-
-  CALL init_dissip
-  CALL init_caldyn
-  CALL init_guided
-  CALL init_advect_tracer
-  CALL init_physics
-!========================================= INITIALIZATION 
-! CALL dynetat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
-  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
-  CALL writefield("phis",f_phis,once=.TRUE.)
-  CALL transfert_request(f_q,req_i1) 
-
+  SUBROUTINE timeloop
+  USE icosa
+  USE dissip_gcm_mod
+  USE caldyn_mod
+  USE etat0_mod
+  USE guided_mod
+  USE advect_tracer_mod
+  USE physics_mod
+  USE mpipara
+  USE omp_para
+  USE trace
+  USE transfert_mod
+  USE check_conserve_mod
+  IMPLICIT NONE  
+    REAL(rstd),POINTER :: phis(:)
+    REAL(rstd),POINTER :: q(:,:,:)
+    REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:)
+    REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:)
+    REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
+    REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:)
+    REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
+    REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
+    REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
+
+    INTEGER :: ind
+    INTEGER :: it,i,j,n,  stage
+    CHARACTER(len=255) :: scheme_name
+    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
+    LOGICAL, PARAMETER :: check=.FALSE.
+
+!$OMP BARRIER
   DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      rhodz=f_rhodz(ind); ps=f_ps(ind)
-     CALL compute_rhodz(.TRUE., ps,rhodz) ! save rhodz for transport scheme before dynamics update ps
+     CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
   END DO
   fluxt_zero=.TRUE.
 
-  ! check that rhodz is consistent with ps
-  CALL transfert_request(f_rhodz,req_i1)
-  CALL transfert_request(f_ps,req_i1)
-  DO ind=1,ndomain
-     CALL swap_dimensions(ind)
-     CALL swap_geometry(ind)
-     rhodz=f_rhodz(ind); ps=f_ps(ind)
-     CALL compute_rhodz(.FALSE., ps, rhodz)   
-  END DO
-
-  CALL transfert_request(f_phis,req_i0) 
-  CALL transfert_request(f_phis,req_i1) 
-  CALL transfert_request(f_phis,req_i1) 
+! check that rhodz is consistent with ps
+!  CALL transfert_request(f_rhodz,req_i1)
+!  CALL transfert_request(f_ps,req_i1)
+!  DO ind=1,ndomain
+!     CALL swap_dimensions(ind)
+!     CALL swap_geometry(ind)
+!     rhodz=f_rhodz(ind); ps=f_ps(ind)
+!     CALL compute_rhodz(.FALSE., ps, rhodz)   
+!  END DO
+  
 
   DO it=0,itaumax
     IF (MOD(it,itau_sync)==0) THEN
-      CALL transfert_request(f_ps,req_i0)
-      CALL transfert_request(f_theta_rhodz,req_i0) 
-      CALL transfert_request(f_u,req_e0_vect)
-      CALL transfert_request(f_q,req_i0) 
+      CALL send_message(f_ps,req_ps0)
+      CALL send_message(f_theta_rhodz,req_theta_rhodz0) 
+      CALL send_message(f_u,req_u0)
+      CALL send_message(f_q,req_q0) 
+      CALL wait_message(req_ps0)
+      CALL wait_message(req_theta_rhodz0) 
+      CALL wait_message(req_u0)
+      CALL wait_message(req_q0) 
     ENDIF
     
+!    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
     IF (mod(it,itau_out)==0 ) THEN
-!      IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
+      CALL writefield("q",f_q)
       CALL update_time_counter(dt*it)
-      CALL compute_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
+      CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
     ENDIF
-
-      CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
+    
+    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
 
     DO stage=1,nb_stage
@@ -246 +252 @@
        CASE (mlf)
           CALL  leapfrog_matsuno_scheme(stage)
+          
           !      CASE ('leapfrog')
-          !      CALL leapfrog_scheme
+          !        CALL leapfrog_scheme
           !
           !      CASE ('adam_bashforth')
-          !      CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
-          !      CALL adam_bashforth_scheme
+          !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
+          !        CALL adam_bashforth_scheme
        CASE DEFAULT
           STOP
@@ -263 +270 @@
 
     IF(MOD(it+1,itau_adv)==0) THEN
-!       CALL transfert_request(f_wfluxt,req_i1) ! FIXME
-!      CALL transfert_request(f_hfluxt,req_e1) ! FIXME
 
        CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
@@ -271 +276 @@
        ! FIXME : check that rhodz is consistent with ps
        IF (check) THEN
-         CALL transfert_request(f_rhodz,req_i1)
-         CALL transfert_request(f_ps,req_i1)
-         CALL transfert_request(f_dps,req_i1) ! FIXME
-         CALL transfert_request(f_wflux,req_i1) ! FIXME
          DO ind=1,ndomain
             CALL swap_dimensions(ind)
             CALL swap_geometry(ind)
-            rhodz=f_rhodz(ind); ps=f_ps(ind); dps=f_dps(ind); 
-            wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
+            rhodz=f_rhodz(ind); ps=f_ps(ind);
             CALL compute_rhodz(.FALSE., ps, rhodz)   
          END DO
        ENDIF
+
     END IF
+
+
+
 !----------------------------------------------------
-  jD_cur = jD_ref + day_ini - day_ref + it/day_step
-  jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step)
-  jD_cur = jD_cur + int(jH_cur)
-  jH_cur = jH_cur - int(jH_cur)
-!       print*,"Just b4 phys"
-    CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
-!----------------------------------------------------
+!    jD_cur = jD_ref + day_ini - day_ref + it/day_step
+!    jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step)
+!    jD_cur = jD_cur + int(jH_cur)
+!    jH_cur = jH_cur - int(jH_cur)
+!    CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
+
 !    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
     ENDDO
+
  
   CONTAINS
@@ -310 +314 @@
          ps=f_ps(ind) ; dps=f_dps(ind) ; 
 
-         DO j=jj_begin,jj_end
-           DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             ps(ij)=ps(ij)+dt*dps(ij)
-           ENDDO
-         ENDDO
+         IF (omp_first) THEN
+           DO j=jj_begin,jj_end
+             DO i=ii_begin,ii_end
+               ij=(j-1)*iim+i
+               ps(ij)=ps(ij)+dt*dps(ij)
+             ENDDO
+           ENDDO
+         ENDIF
+         
          hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
@@ -324 +331 @@
        du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
 
-       DO l=1,llm
+       DO l=ll_begin,ll_end
          DO j=jj_begin,jj_end
            DO i=ii_begin,ii_end
@@ -342 +349 @@
 
     SUBROUTINE RK_scheme(stage)
+      USE caldyn_gcm_mod
       IMPLICIT NONE
       INTEGER :: ind, stage
@@ -351 +359 @@
 
       tau = dt*coef(stage)
+
+      DO ind=1,ndomain
+         CALL swap_dimensions(ind)
+         CALL swap_geometry(ind)
+         ps=f_ps(ind)   
+         psm1=f_psm1(ind) 
+         dps=f_dps(ind) 
+         
+         IF (stage==1) THEN ! first stage : save model state in XXm1
+           IF (omp_first) THEN
+             DO j=jj_begin,jj_end
+               DO i=ii_begin,ii_end
+                 ij=(j-1)*iim+i
+                 psm1(ij)=ps(ij)
+               ENDDO
+             ENDDO
+           ENDIF
+         END IF
+         
+         ! updates are of the form x1 := x0 + tau*f(x1)
+           IF (omp_first) THEN
+             DO j=jj_begin,jj_end
+               DO i=ii_begin,ii_end
+                 ij=(j-1)*iim+i
+                 ps(ij)=psm1(ij)+tau*dps(ij)
+               ENDDO
+             ENDDO
+           ENDIF
+       
+       ENDDO
+
+       CALL send_message(f_ps,req_ps)      
+
+
       
       DO ind=1,ndomain
@@ -360 +402 @@
          
          IF (stage==1) THEN ! first stage : save model state in XXm1
-
-           DO j=jj_begin,jj_end
-             DO i=ii_begin,ii_end
-               ij=(j-1)*iim+i
-               psm1(ij)=ps(ij)
-             ENDDO
-           ENDDO
               
-           DO l=1,llm
+         DO l=ll_begin,ll_end
              DO j=jj_begin,jj_end
                DO i=ii_begin,ii_end
@@ -382 +417 @@
          END IF
          ! updates are of the form x1 := x0 + tau*f(x1)
-         DO j=jj_begin,jj_end
-           DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             ps(ij)=psm1(ij)+tau*dps(ij)
-           ENDDO
-         ENDDO
          
-         DO l=1,llm
+         DO l=ll_begin,ll_end
            DO j=jj_begin,jj_end
              DO i=ii_begin,ii_end
@@ -468 +497 @@
 
         
-!        IF (MOD(it,matsuno_period+1)==0) THEN
         IF (stage==1) THEN
           psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
@@ -475 +503 @@
           theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
 
-!        ELSE IF (MOD(it,matsuno_period+1)==1) THEN
         ELSE IF (stage==2) THEN
 
@@ -538 +565 @@
     USE icosa
     USE disvert_mod
+    USE omp_para
     LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
     REAL(rstd), INTENT(IN) :: ps(iim*jjm)
@@ -544 +572 @@
     INTEGER :: l,i,j,ij,dd
     err=0.
-    IF(comp) THEN 
-       dd=1
-    ELSE
-!       dd=-1
-       dd=0
-    END IF
-
-    DO l = 1, llm
+
+    dd=0
+
+    DO l = ll_begin, ll_end
        DO j=jj_begin-dd,jj_end+dd
           DO i=ii_begin-dd,ii_end+dd
              ij=(j-1)*iim+i
-             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
+             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
              IF(comp) THEN
                 rhodz(ij,l) = m
@@ -570 +594 @@
           STOP
        ELSE
-          PRINT *, 'No discrepancy between ps and rhodz detected'
+!          PRINT *, 'No discrepancy between ps and rhodz detected'
        END IF
     END IF
@@ -578 +602 @@
   SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
     USE icosa
+    USE omp_para
     REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
     REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1)
@@ -585 +610 @@
 
     IF(fluxt_zero) THEN
-!       PRINT *, 'Accumulating fluxes (first)'
+
        fluxt_zero=.FALSE.
-       DO l=1,llm
+
+       DO l=ll_begin,ll_end
          DO j=jj_begin-1,jj_end+1
            DO i=ii_begin-1,ii_end+1
@@ -598 +624 @@
        ENDDO
 
-       DO l=1,llm+1
+       DO l=ll_begin,ll_endp1
          DO j=jj_begin,jj_end
            DO i=ii_begin,ii_end
@@ -607 +633 @@
        ENDDO
     ELSE
-!       PRINT *, 'Accumulating fluxes (next)'
-       DO l=1,llm
+
+       DO l=ll_begin,ll_end
          DO j=jj_begin-1,jj_end+1
            DO i=ii_begin-1,ii_end+1
@@ -619 +645 @@
        ENDDO
 
-       DO l=1,llm+1
+       DO l=ll_begin,ll_endp1
          DO j=jj_begin,jj_end
            DO i=ii_begin,ii_end
@@ -629 +655 @@
 
     END IF
+
   END SUBROUTINE accumulate_fluxes
   

codes/icosagcm/trunk/src/trace.F90

@@ -1 +1 @@
 MODULE trace
 
-
+  INTEGER,SAVE :: markId
+  
 CONTAINS
 
+  SUBROUTINE init_trace
+  IMPLICIT NONE
+#ifdef VTRACE
+#include <vt_user.inc>
+#endif 
+
+#ifdef VTRACE
+     VT_MARKER_DEF("marker", VT_MARKER_TYPE_HINT, markId)
+#endif
+  
+  END SUBROUTINE init_trace
+  
+  
   SUBROUTINE trace_start(name)
   IMPLICIT NONE
@@ -11 +25 @@
 #endif 
 
+!$OMP MASTER
 #ifdef VTRACE
      VT_USER_START(name)
 #endif
+!$OMP END MASTER
 
   END SUBROUTINE trace_start    
@@ -25 +41 @@
     CHARACTER(LEN=*),INTENT(IN) :: name
 
+!$OMP MASTER
 #ifdef VTRACE
      VT_USER_END(name)
 #endif
+!$OMP END MASTER
 
   END SUBROUTINE trace_end    
 
+  SUBROUTINE Marker(name)
+  IMPLICIT NONE
+  CHARACTER(LEN=*),INTENT(IN) :: name
+#ifdef VTRACE
+#include <vt_user.inc>
+#endif 
+
+!$OMP MASTER
+#ifdef VTRACE
+     VT_MARKER(markId,name)
+#endif
+!$OMP END MASTER
+
+  END SUBROUTINE Marker
+
 END MODULE trace

codes/icosagcm/trunk/src/transfert.F90

@@ -3 +3 @@
 #ifdef CPP_USING_MPI
   USE transfert_mpi_mod, ONLY : init_transfert, transfert_request=>transfert_request_mpi, req_i1,req_e1_vect, &
-                                req_e1_scal, req_i0, req_e0_vect, req_e0_scal, request_add_point, create_request, gather_field
+                                req_e1_scal, req_i0, req_e0_vect, req_e0_scal, request_add_point, create_request, gather_field, &
+                                t_message,init_message=>init_message_mpi,transfert_message=>transfert_message_mpi,  &
+                                send_message=>send_message_mpi,test_message=>test_message_mpi,wait_message=>wait_message_mpi,barrier
 #else 
-  USE transfert_mpi_mod, ONLY : init_transfert, transfert_request, req_i1,req_e1_vect, &
-                                req_e1_scal,request_add_point, create_request, gather_field
+  USE transfert_mpi_mod, ONLY : init_transfert, transfert_request=>transfert_request_seq, req_i1,req_e1_vect, &
+                                req_e1_scal,req_i0, req_e0_vect, req_e0_scal, request_add_point, create_request, gather_field, &
+                                t_message,init_message=>init_message_seq,transfert_message=>transfert_message_seq,  &
+                                send_message=>send_message_seq,test_message=>test_message_seq,wait_message=>wait_message_seq,barrier
 #endif
   

codes/icosagcm/trunk/src/transfert_mpi.f90

@@ -1 +1 @@
 MODULE transfert_mpi_mod
 USE genmod
+USE field_mod
   
   TYPE array
@@ -13 +14 @@
     REAL,POINTER    :: buffer_r4(:,:,:)
   END TYPE array
+  
+  TYPE t_buffer
+    REAL,POINTER    :: r2(:)
+    REAL,POINTER    :: r3(:,:)
+    REAL,POINTER    :: r4(:,:,:)
+  END TYPE t_buffer    
     
   TYPE t_request
@@ -41 +48 @@
   TYPE(t_request),POINTER :: req_e0_vect(:)
   
+  TYPE t_message
+    TYPE(t_request), POINTER :: request(:)
+    INTEGER :: nreq
+    INTEGER, POINTER :: mpi_req(:)
+    INTEGER, POINTER :: status(:,:)
+    TYPE(t_buffer),POINTER :: buffers(:) 
+    TYPE(t_field),POINTER :: field(:)
+    LOGICAL :: completed
+    LOGICAL :: pending
+  END TYPE t_message
   
 CONTAINS
-
+  
   SUBROUTINE init_transfert
   USE domain_mod
@@ -70 +87 @@
       DO j=jj_begin,jj_end+1
         CALL request_add_point(ind,ii_begin-1,j,req_i1)
-      ENDDO    
-    
-      DO i=ii_begin,ii_end
-!        CALL request_add_point(ind,i,jj_begin,req_i1)
-      ENDDO
-
-      DO j=jj_begin,jj_end
-!        CALL request_add_point(ind,ii_end,j,req_i1)
-      ENDDO    
-    
-      DO i=ii_begin,ii_end
-!        CALL request_add_point(ind,i,jj_end,req_i1)
-      ENDDO    
-
-      DO j=jj_begin,jj_end
-!        CALL request_add_point(ind,ii_begin,j,req_i1)
       ENDDO    
     
@@ -142 +143 @@
       ENDDO   
 
-      DO i=ii_begin+1,ii_end-1
-!        CALL request_add_point(ind,i,jj_begin,req_e1_scal,right)
-!        CALL request_add_point(ind,i,jj_end,req_e1_scal,right)
-      ENDDO
-    
-      DO j=jj_begin+1,jj_end-1
-!        CALL request_add_point(ind,ii_begin,j,req_e1_scal,rup)
-!        CALL request_add_point(ind,ii_end,j,req_e1_scal,rup)
-      ENDDO   
-
-!      CALL request_add_point(ind,ii_begin+1,jj_begin,req_e1_scal,left)
-!      CALL request_add_point(ind,ii_begin,jj_begin+1,req_e1_scal,ldown)
-!      CALL request_add_point(ind,ii_begin+1,jj_end,req_e1_scal,left)
-!      CALL request_add_point(ind,ii_end,jj_begin+1,req_e1_scal,ldown)
-
     ENDDO
 
@@ -211 +197 @@
       ENDDO   
 
-      DO i=ii_begin+1,ii_end-1
-!        CALL request_add_point(ind,i,jj_begin,req_e1_vect,right)
-!        CALL request_add_point(ind,i,jj_end,req_e1_vect,right)
-      ENDDO
-    
-      DO j=jj_begin+1,jj_end-1
-!        CALL request_add_point(ind,ii_begin,j,req_e1_vect,rup)
-!        CALL request_add_point(ind,ii_end,j,req_e1_vect,rup)
-      ENDDO   
-
-!      CALL request_add_point(ind,ii_begin+1,jj_begin,req_e1_vect,left)
-!      CALL request_add_point(ind,ii_begin,jj_begin+1,req_e1_vect,ldown)
-!      CALL request_add_point(ind,ii_begin+1,jj_end,req_e1_vect,left)
-!      CALL request_add_point(ind,ii_end,jj_begin+1,req_e1_vect,ldown)
   
     ENDDO  
@@ -322 +294 @@
     target_j=>req%target_j
     target_sign=>req%target_sign
-!    req%max_size=req%max_size*2
+
     ALLOCATE(req%src_domain(req%max_size*2))
     ALLOCATE(req%src_ind(req%max_size*2))
@@ -610 +582 @@
 
 
+  SUBROUTINE init_message_seq(field, request, message)
+  USE field_mod
+  USE domain_mod
+  USE mpi_mod
+  USE mpipara
+  USE mpi_mod
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: field(:)
+    TYPE(t_request),POINTER :: request(:)
+    TYPE(t_message) :: message
+
+!$OMP MASTER    
+    message%request=>request
+!$OMP END MASTER    
+!$OMP BARRIER    
+
+  END SUBROUTINE init_message_seq
+
+  SUBROUTINE send_message_seq(field,message)
+  USE field_mod
+  USE domain_mod
+  USE mpi_mod
+  USE mpipara
+  USE omp_para
+  USE trace
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: field(:)
+    TYPE(t_message) :: message
+
+    CALL transfert_request_seq(field,message%request)
+    
+  END SUBROUTINE send_message_seq
+  
+  SUBROUTINE test_message_seq(message)
+  IMPLICIT NONE
+    TYPE(t_message) :: message
+  END SUBROUTINE  test_message_seq
+  
+   
+  SUBROUTINE wait_message_seq(message)
+  IMPLICIT NONE
+    TYPE(t_message) :: message
+    
+  END SUBROUTINE wait_message_seq    
+
+  SUBROUTINE transfert_message_seq(field,message)
+  USE field_mod
+  USE domain_mod
+  USE mpi_mod
+  USE mpipara
+  USE omp_para
+  USE trace
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: field(:)
+    TYPE(t_message) :: message
+
+   CALL send_message_seq(field,message)
+    
+  END SUBROUTINE transfert_message_seq    
+    
+    
+  SUBROUTINE init_message_mpi(field,request, message)
+  USE field_mod
+  USE domain_mod
+  USE mpi_mod
+  USE mpipara
+  USE mpi_mod
+  IMPLICIT NONE
+  
+    TYPE(t_field),POINTER :: field(:)
+    TYPE(t_request),POINTER :: request(:)
+    TYPE(t_message) :: message
+
+    TYPE(ARRAY),POINTER :: recv,send 
+    TYPE(t_request),POINTER :: req
+    INTEGER :: irecv,isend
+    INTEGER :: ireq,nreq
+    INTEGER :: ind
+    INTEGER :: dim3,dim4
+
+!$OMP MASTER
+    message%request=>request
+    nreq=sum(request(:)%nsend)+sum(request(:)%nrecv)
+    message%nreq=nreq
+    ALLOCATE(message%mpi_req(nreq))
+    ALLOCATE(message%buffers(nreq))
+    ALLOCATE(message%status(MPI_STATUS_SIZE,nreq))
+    
+    message%pending=.FALSE.
+    message%completed=.FALSE.
+    
+    IF (field(1)%data_type==type_real) THEN
+
+      IF (field(1)%ndim==2) THEN
+      
+        ireq=0
+        DO ind=1,ndomain
+          req=>request(ind)
+      
+          DO isend=1,req%nsend
+            ireq=ireq+1
+            send=>req%send(isend)
+            CALL allocate_mpi_buffer(message%buffers(ireq)%r2,send%size)
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            CALL allocate_mpi_buffer(message%buffers(ireq)%r2,recv%size)
+          ENDDO
+        
+        ENDDO
+      
+      
+      ELSE  IF (field(1)%ndim==3) THEN
+    
+        ireq=0
+        DO ind=1,ndomain
+          dim3=size(field(ind)%rval3d,2)
+          req=>request(ind)
+ 
+          DO isend=1,req%nsend
+            ireq=ireq+1
+            send=>req%send(isend)
+            CALL allocate_mpi_buffer(message%buffers(ireq)%r3,send%size,dim3)
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            CALL allocate_mpi_buffer(message%buffers(ireq)%r3,recv%size,dim3)
+
+          ENDDO
+        
+        ENDDO
+
+
+      ELSE  IF (field(1)%ndim==4) THEN
+   
+        ireq=0
+        DO ind=1,ndomain
+          dim3=size(field(ind)%rval4d,2)
+          dim4=size(field(ind)%rval4d,3)
+          req=>request(ind)
+
+          DO isend=1,req%nsend
+            ireq=ireq+1
+            send=>req%send(isend)
+            CALL allocate_mpi_buffer(message%buffers(ireq)%r4,send%size,dim3,dim4)
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            CALL allocate_mpi_buffer(message%buffers(ireq)%r4,recv%size,dim3,dim4)
+          ENDDO
+        
+        ENDDO
+      
+      ENDIF      
+    ENDIF
+!$OMP END MASTER
+!$OMP BARRIER    
+  END SUBROUTINE init_message_mpi
+  
+  SUBROUTINE barrier
+  USE mpi_mod
+  USE mpipara
+  IMPLICIT NONE
+    
+    CALL MPI_BARRIER(comm_icosa,ierr)
+    
+  END SUBROUTINE barrier  
+    
+  SUBROUTINE transfert_message_mpi(field,message)
+  USE field_mod
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: field(:)
+    TYPE(t_message) :: message
+    
+    CALL send_message_mpi(field,message)
+    CALL wait_message_mpi(message)
+    
+  END SUBROUTINE transfert_message_mpi
+
+
+  SUBROUTINE send_message_mpi(field,message)
+  USE field_mod
+  USE domain_mod
+  USE mpi_mod
+  USE mpipara
+  USE omp_para
+  USE trace
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: field(:)
+    TYPE(t_message) :: message
+    REAL(rstd),POINTER :: rval2d(:) 
+    REAL(rstd),POINTER :: rval3d(:,:) 
+    REAL(rstd),POINTER :: rval4d(:,:,:) 
+    REAL(rstd),POINTER :: buffer_r2(:) 
+    REAL(rstd),POINTER :: buffer_r3(:,:) 
+    REAL(rstd),POINTER :: buffer_r4(:,:,:) 
+    INTEGER,POINTER :: value(:) 
+    INTEGER,POINTER :: sgn(:) 
+    TYPE(ARRAY),POINTER :: recv,send 
+    TYPE(t_request),POINTER :: req
+    INTEGER, ALLOCATABLE :: mpi_req(:)
+    INTEGER, ALLOCATABLE :: status(:,:)
+    INTEGER :: irecv,isend
+    INTEGER :: ireq,nreq
+    INTEGER :: ind,n,l,m
+    INTEGER :: dim3,dim4
+
+!$OMP BARRIER
+
+    CALL trace_start("transfert_mpi")
+
+    nreq=message%nreq
+    message%field=>field
+
+!$OMP MASTER
+    message%completed=.FALSE.
+    message%pending=.TRUE.
+!$OMP END MASTER
+    
+    IF (field(1)%data_type==type_real) THEN
+      IF (field(1)%ndim==2) THEN
+
+        ireq=0
+        DO ind=1,ndomain
+          rval2d=>field(ind)%rval2d
+        
+          req=>message%request(ind)
+          DO isend=1,req%nsend
+            ireq=ireq+1
+            send=>req%send(isend)
+            buffer_r2=>message%buffers(ireq)%r2
+            value=>send%value
+
+            CALL trace_in
+
+!$OMP DO SCHEDULE(STATIC)
+            DO n=1,send%size
+              buffer_r2(n)=rval2d(value(n))
+            ENDDO
+            
+            CALL trace_out
+
+!$OMP MASTER
+            CALL MPI_ISSEND(buffer_r2,send%size,MPI_REAL8,send%rank,ind,comm_icosa, message%mpi_req(ireq),ierr)
+!$OMP END MASTER
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            buffer_r2=>message%buffers(ireq)%r2
+!$OMP MASTER
+            CALL MPI_IRECV(buffer_r2,recv%size,MPI_REAL8,recv%rank,recv%domain,comm_icosa, message%mpi_req(ireq),ierr)
+!$OMP END MASTER
+          ENDDO
+        
+        ENDDO
+      
+      ELSE  IF (field(1)%ndim==3) THEN
+      
+        ireq=0
+        DO ind=1,ndomain
+          dim3=size(field(ind)%rval3d,2)
+          rval3d=>field(ind)%rval3d
+          req=>message%request(ind)
+ 
+          DO isend=1,req%nsend
+            ireq=ireq+1
+            send=>req%send(isend)
+            buffer_r3=>message%buffers(ireq)%r3
+            value=>send%value
+
+            CALL trace_in
+            
+!$OMP DO SCHEDULE(STATIC)
+              DO n=1,send%size
+                buffer_r3(n,:)=rval3d(value(n),:)
+              ENDDO
+
+             CALL trace_out
+
+!$OMP MASTER
+            CALL MPI_ISSEND(buffer_r3,send%size*dim3,MPI_REAL8,send%rank,ind,comm_icosa, message%mpi_req(ireq),ierr)
+!$OMP END MASTER
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            buffer_r3=>message%buffers(ireq)%r3
+!$OMP MASTER           
+            CALL MPI_IRECV(buffer_r3,recv%size*dim3,MPI_REAL8,recv%rank,recv%domain,comm_icosa, message%mpi_req(ireq),ierr)
+!$OMP END MASTER
+
+          ENDDO
+        
+        ENDDO
+
+      ELSE  IF (field(1)%ndim==4) THEN
+    
+        ireq=0
+        DO ind=1,ndomain
+          dim3=size(field(ind)%rval4d,2)
+          dim4=size(field(ind)%rval4d,3)
+          rval4d=>field(ind)%rval4d
+          req=>message%request(ind)
+
+          DO isend=1,req%nsend
+            ireq=ireq+1
+            send=>req%send(isend)
+            buffer_r4=>message%buffers(ireq)%r4
+            value=>send%value
+
+            CALL trace_in
+
+!$OMP DO SCHEDULE(STATIC)
+            DO n=1,send%size
+               buffer_r4(n,:,:)=rval4d(value(n),:,:)
+            ENDDO
+
+           CALL trace_out
+
+!$OMP MASTER
+            CALL MPI_ISSEND(buffer_r4,send%size*dim3*dim4,MPI_REAL8,send%rank,ind,comm_icosa, message%mpi_req(ireq),ierr)
+!$OMP END MASTER
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            buffer_r4=>message%buffers(ireq)%r4
+!$OMP MASTER           
+            CALL MPI_IRECV(buffer_r4,recv%size*dim3*dim4,MPI_REAL8,recv%rank,recv%domain,comm_icosa, message%mpi_req(ireq),ierr)
+!$OMP END MASTER
+          ENDDO
+        
+        ENDDO
+      
+      ENDIF      
+      
+    ENDIF
+
+    CALL trace_end("transfert_mpi")
+!$OMP BARRIER
+    
+  END SUBROUTINE send_message_mpi
+  
+  SUBROUTINE test_message_mpi(message)
+  IMPLICIT NONE
+    TYPE(t_message) :: message
+    
+    INTEGER :: ierr
+!$OMP MASTER
+     IF (.NOT. message%pending) RETURN
+!$OMP END MASTER
+
+!$OMP MASTER
+     IF (.NOT. message%completed) CALL MPI_TESTALL(message%nreq,message%mpi_req,message%completed,message%status,ierr)
+!$OMP END MASTER
+  END SUBROUTINE  test_message_mpi
+  
+   
+  SUBROUTINE wait_message_mpi(message)
+  USE field_mod
+  USE domain_mod
+  USE mpi_mod
+  USE mpipara
+  USE omp_para
+  USE trace
+  IMPLICIT NONE
+    TYPE(t_message) :: message
+
+    TYPE(t_field),POINTER :: field(:)
+    REAL(rstd),POINTER :: rval2d(:) 
+    REAL(rstd),POINTER :: rval3d(:,:) 
+    REAL(rstd),POINTER :: rval4d(:,:,:) 
+    REAL(rstd),POINTER :: buffer_r2(:) 
+    REAL(rstd),POINTER :: buffer_r3(:,:) 
+    REAL(rstd),POINTER :: buffer_r4(:,:,:) 
+    INTEGER,POINTER :: value(:) 
+    INTEGER,POINTER :: sgn(:) 
+    TYPE(ARRAY),POINTER :: recv,send 
+    TYPE(t_request),POINTER :: req
+    INTEGER, ALLOCATABLE :: mpi_req(:)
+    INTEGER, ALLOCATABLE :: status(:,:)
+    INTEGER :: irecv,isend
+    INTEGER :: ireq,nreq
+    INTEGER :: ind,n,l,m
+    INTEGER :: dim3,dim4
+
+!$OMP BARRIER
+
+    CALL trace_start("transfert_mpi")
+
+    IF (.NOT. message%pending) RETURN
+    
+    field=>message%field
+    nreq=message%nreq
+    
+    IF (field(1)%data_type==type_real) THEN
+      IF (field(1)%ndim==2) THEN
+
+!$OMP MASTER
+        IF (.NOT. message%completed) CALL MPI_WAITALL(nreq,message%mpi_req,message%status,ierr)
+!$OMP END MASTER
+!$OMP BARRIER
+
+        ireq=0        
+        DO ind=1,ndomain
+          rval2d=>field(ind)%rval2d
+          req=>message%request(ind)
+
+          DO isend=1,req%nsend
+            ireq=ireq+1
+          ENDDO
+     
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            buffer_r2=>message%buffers(ireq)%r2
+            value=>recv%value
+            sgn=>recv%sign
+
+            CALL trace_in
+            
+!$OMP DO SCHEDULE(STATIC)
+            DO n=1,recv%size
+              rval2d(value(n))=buffer_r2(n)*sgn(n)  
+            ENDDO        
+
+            CALL trace_out
+
+          ENDDO
+        
+        ENDDO
+      
+      
+      ELSE  IF (field(1)%ndim==3) THEN
+
+!$OMP MASTER
+        IF (.NOT. message%completed) CALL MPI_WAITALL(nreq,message%mpi_req,message%status,ierr)
+!$OMP END MASTER
+!$OMP BARRIER
+
+        ireq=0        
+        DO ind=1,ndomain
+          rval3d=>field(ind)%rval3d
+          req=>message%request(ind)
+
+          DO isend=1,req%nsend
+            ireq=ireq+1
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            buffer_r3=>message%buffers(ireq)%r3
+            value=>recv%value
+            sgn=>recv%sign
+
+            CALL trace_in
+            
+!$OMP DO SCHEDULE(STATIC)
+            DO n=1,recv%size
+              rval3d(value(n),:)=buffer_r3(n,:)*sgn(n)  
+            ENDDO  
+
+            CALL trace_out
+
+          ENDDO
+        
+        ENDDO
+
+      ELSE  IF (field(1)%ndim==4) THEN
+!$OMP MASTER
+        IF (.NOT. message%completed) CALL MPI_WAITALL(nreq,message%mpi_req,message%status,ierr)
+!$OMP END MASTER
+!$OMP BARRIER
+
+        ireq=0        
+        DO ind=1,ndomain
+          rval4d=>field(ind)%rval4d
+          req=>message%request(ind)
+
+          DO isend=1,req%nsend
+            ireq=ireq+1
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            buffer_r4=>message%buffers(ireq)%r4
+            value=>recv%value
+            sgn=>recv%sign
+
+            CALL trace_in
+
+!$OMP DO SCHEDULE(STATIC)
+            DO n=1,recv%size
+              rval4d(value(n),:,:)=buffer_r4(n,:,:)*sgn(n) 
+            ENDDO
+
+            CALL trace_out
+
+          ENDDO
+        
+        ENDDO
+      
+      ENDIF      
+      
+    ENDIF
+
+!$OMP MASTER
+    message%pending=.FALSE.
+!$OMP END MASTER
+
+    CALL trace_end("transfert_mpi")
+!$OMP BARRIER
+    
+  END SUBROUTINE wait_message_mpi
+
+
   SUBROUTINE transfert_request_mpi(field,request)
   USE field_mod
@@ -637 +1137 @@
 
     CALL trace_start("transfert_mpi")
-    
+
     IF (field(1)%data_type==type_real) THEN
       IF (field(1)%ndim==2) THEN
@@ -675 +1175 @@
         
         CALL MPI_WAITALL(nreq,mpi_req,status,ierr)
-!        DO ind=1,ndomain
-!          field(ind)%rval2d(:)=0.
-!        ENDDO
         
         DO ind=1,ndomain
@@ -739 +1236 @@
         
         CALL MPI_WAITALL(nreq,mpi_req,status,ierr)
-!        DO ind=1,ndomain
-!          field(ind)%rval2d(:)=0.
-!        ENDDO
         
         DO ind=1,ndomain
@@ -803 +1297 @@
         
         CALL MPI_WAITALL(nreq,mpi_req,status,ierr)
-!        DO ind=1,ndomain
-!          field(ind)%rval2d(:)=0.
-!        ENDDO
         
         DO ind=1,ndomain
@@ -837 +1328 @@
   END SUBROUTINE transfert_request_mpi
    
-  SUBROUTINE transfert_request(field,request)
+  SUBROUTINE transfert_request_seq(field,request)
   USE field_mod
   USE domain_mod
@@ -875 +1366 @@
     ENDDO
     
-  END SUBROUTINE transfert_request
+  END SUBROUTINE transfert_request_seq
   
   
@@ -949 +1440 @@
         
   END SUBROUTINE gather_field
-  
+
+   
+  SUBROUTINE trace_in
+  USE trace
+  IMPLICIT NONE
+  
+    CALL trace_start("transfert_buffer")
+  END SUBROUTINE trace_in              
+
+  SUBROUTINE trace_out
+  USE trace
+  IMPLICIT NONE
+  
+    CALL trace_end("transfert_buffer")
+  END SUBROUTINE trace_out              
 
 END MODULE transfert_mpi_mod

codes/icosagcm/trunk/src/wind.f90

@@ -3 +3 @@
 CONTAINS
 
+  SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat)
+  USE icosa
+  IMPLICIT NONE
+    TYPE(t_field), POINTER :: f_u(:) ! IN  : normal velocity components on edges
+    TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons
+    
+    REAL(rstd),POINTER :: u(:,:),  ulon(:,:), ulat(:,:)
+    INTEGER :: ind
+
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       u=f_u(ind)
+       ulon=f_ulon(ind)
+       ulat=f_ulat(ind)
+       CALL compute_un2ulonlat(u,ulon, ulat)
+    END DO
+
+  END SUBROUTINE un2ulonlat
 
  
@@ -214 +233 @@
  END SUBROUTINE compute_wind_centered_lonlat_compound
 
+ SUBROUTINE compute_un2ulonlat(un, ulon, ulat)
+  USE icosa  
+    
+  IMPLICIT NONE
+  REAL(rstd),INTENT(IN)  :: un(3*iim*jjm,llm)
+  REAL(rstd),INTENT(OUT) :: ulon(iim*jjm,llm)
+  REAL(rstd),INTENT(OUT) :: ulat(iim*jjm,llm)
+
+  REAL(rstd)             :: uc(iim*jjm,3,llm)
+    
+  CALL compute_wind_centered(un,uc)  
+  CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
+ 
+ END SUBROUTINE compute_un2ulonlat
 
 END MODULE wind_mod

codes/icosagcm/trunk/src/write_field.f90

@@ -75 +75 @@
       
       TYPE(t_field),POINTER :: field_glo(:)
-      
+
+!$OMP MASTER      
       IF(PRESENT(once)) THEN
          once_=once
@@ -101 +102 @@
       
       CALL deallocate_field_glo(field_glo)
+!$OMP END MASTER
       
    END SUBROUTINE writefield
@@ -333 +335 @@
     USE domain_mod
     USE field_mod
-!    USE dimensions
-!    USE geometry
     IMPLICIT NONE
       CHARACTER(LEN=*),INTENT(IN) :: name_in
@@ -482 +482 @@
         ENDIF
         
-!        PRINT *,NF90_STRERROR(status)
-!        ncell=ncell+n
-
-!     ENDDO
      
      ELSE IF (Field(ind_b)%field_type==field_Z) THEN
@@ -604 +600 @@
           ENDDO
         ENDIF
-        
-!        PRINT *,NF90_STRERROR(status)
-!        ncell=ncell+n
-!
-!     ENDDO
      
      ENDIF
@@ -712 +703 @@
             ENDDO
           ENDDO
-!           DO l=1,size(field(ind)%rval3d,2)
-            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,  &
+          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,  &
                                 start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /))
-!           ENDDO
           DEALLOCATE(field_val3d)
+
         ELSE IF (field(1)%ndim==4) THEN
 
@@ -732 +722 @@
               ENDDO
              ENDDO
-!            DO l=1,size(field(ind)%rval4d,2)
+
             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d(:,l), &
                                 start=(/ displ+ncell,l,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /))
-!            ENDDO
             DEALLOCATE(field_val3d)
           ENDDO         
         ENDIF
         
-!        PRINT *,NF90_STRERROR(status)
         ncell=ncell+n
 
@@ -810 +798 @@
             ENDDO
           ENDDO
-!           DO l=1,size(field(ind)%rval3d,2)
+
            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,    &
                                start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /))
-!           ENDDO
           DEALLOCATE(field_val3d)
+
         ELSE IF (field(1)%ndim==4) THEN
 
@@ -836 +824 @@
             ENDDO
 
-!            DO l=1,size(field(ind)%rval4d,2)
-
             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,  &
                                 start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /))
-!            ENDDO
             DEALLOCATE(field_val3d)
           ENDDO
         ENDIF
         
-!        PRINT *,NF90_STRERROR(status)
         ncell=ncell+n
 
@@ -1159 +1143 @@
     USE metric
     USE spherical_geom_mod
-!    USE dimensions
-!    USE geometry
     IMPLICIT NONE
       CHARACTER(LEN=*),INTENT(IN) :: name_in
@@ -1292 +1274 @@
       
         status = NF90_ENDDEF(ncid)      
-
-!        ncell=1
-!        DO ind=ind_b,ind_e
-!          d=>domain_type(ind)
-  
-!          n=0
-!          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-!            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-!              IF (single .OR. d%assign_domain(i,j)==ind) n=n+1
-!            ENDDO
-!          ENDDO
-          
+         
          ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell))
           
@@ -1330 +1301 @@
          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
   
-!         ncell=ncell+n
          DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
-!      END DO 
 
     ELSE IF (Field(ind_b)%field_type==field_Z) THEN
@@ -1416 +1385 @@
         status = NF90_ENDDEF(ncid)      
 
-!        ncell=1
-!        DO ind=ind_b,ind_e
-!          d=>domain_type(ind)
-!  
-!          n=0
-!          DO j=d%jj_begin+1,d%jj_end
-!            DO i=d%ii_begin,d%ii_end-1
-!              n=n+1
-!            ENDDO
-!          ENDDO
-!
-!          DO j=d%jj_begin,d%jj_end-1
-!            DO i=d%ii_begin+1,d%ii_end
-!              n=n+1
-!            ENDDO
-!          ENDDO
-
-!         ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
          ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell))
           
@@ -1446 +1397 @@
                CALL xyz2lonlat(d%vertex(:,vdown,i,j),lon(n),lat(n))
                lon(n)=lon(n)*180/Pi
- !              IF (lon(n)<0) lon(n)=lon(n)+360
                lat(n)=lat(n)*180/Pi
-!               CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
-!               CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
-!               CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
 
                CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n))
@@ -1459 +1406 @@
                  bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
                  bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
                ENDDO
              ENDDO
@@ -1468 +1414 @@
                nij=(j-1)*d%iim+i
                n=n+1
-!              CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
                CALL xyz2lonlat(d%vertex(:,vup,i,j),lon(n),lat(n))
                lon(n)=lon(n)*180/Pi
-!              IF (lon(n)<0) lon(n)=lon(n)+360
                lat(n)=lat(n)*180/Pi
-!              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
-!              CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
-!              CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
                CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n))
                CALL xyz2lonlat(d%xyz(:,i,j+1),bounds_lon(1,n), bounds_lat(1,n))
@@ -1483 +1424 @@
                  bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
                  bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-!                 IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
                ENDDO
              ENDDO
@@ -1494 +1434 @@
          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
   
-!          ncell=ncell+n
           DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
-!      END DO          
-    ENDIF
-
-
-      
-!      status = NF90_CLOSE(ncid)
+      ENDIF
+
 
     END SUBROUTINE Create_Header_gen
@@ -1783 +1718 @@
               CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n))
               lon(n)=lon(n)*180/Pi
- !             IF (lon(n)<0) lon(n)=lon(n)+360
               lat(n)=lat(n)*180/Pi
               CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
@@ -1792 +1726 @@
                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
               ENDDO
             ENDDO
@@ -1803 +1736 @@
               CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
               lon(n)=lon(n)*180/Pi
-!              IF (lon(n)<0) lon(n)=lon(n)+360
               lat(n)=lat(n)*180/Pi
               CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
@@ -1812 +1744 @@
                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
               ENDDO
             ENDDO
@@ -1829 +1760 @@
 
 
-      
-!      status = NF90_CLOSE(ncid)
-
-    end subroutine Create_Header_mpi  
+   END SUBROUTINE Create_Header_mpi  
    
    SUBROUTINE Close_files