Changeset 581

Timestamp:

10/13/17 16:00:58 (7 years ago)

Author:

dubos

Message:

trunk : upgrading to devel

Location:

codes/icosagcm/trunk/src

Files:

• diagnostics/observable.f90 (modified) (9 diffs)
• initial/etat0.f90 (modified) (13 diffs)
• output/write_etat0.f90 (modified) (3 diffs)
• time/timeloop_gcm.f90 (modified) (2 diffs)

codes/icosagcm/trunk/src/diagnostics/observable.f90

@@ -5 +5 @@
 
   TYPE(t_field),POINTER, SAVE :: f_buf_i(:), &
-       f_buf_uh(:), & ! horizontal velocity, different from prognostic velocity if NH
+       f_buf_Fel(:), f_buf_uh(:), & ! horizontal velocity, different from prognostic velocity if NH
        f_buf_ulon(:), f_buf_ulat(:), &
        f_buf_u3d(:) ! unused, remove ?
@@ -28 +28 @@
     CALL allocate_field(f_buf_ulat,field_t,type_real,llm, name="buf_ulat")
     CALL allocate_field(f_buf_uh,  field_u,type_real,llm, name="buf_uh")
+    CALL allocate_field(f_buf_Fel, field_u,type_real,llm+1, name="buf_F_el")
     CALL allocate_field(f_buf_v,   field_z,type_real,llm, name="buf_v")
     CALL allocate_field(f_buf_s,   field_t,type_real, name="buf_s")
@@ -60 +61 @@
     
     IF(init) THEN
+       IF(is_master) PRINT *, 'Creating output files ...'
        scalar(1)=dt
        IF (is_omp_master) CALL xios_send_field("timestep", scalar)
@@ -74 +76 @@
 
        CALL output_field("phis",f_phis)
-       CALL output_field("Ai",geom%Ai)       
+       CALL output_field("Ai",geom%Ai) 
+       IF(is_master) PRINT *, '... done creating output files. Writing initial condition ...'
     END IF
 
@@ -120 +123 @@
        CALL output_field("geopot_init",f_geopot)
        CALL output_field("q_init",f_q)
+       IF(is_master) PRINT *, 'Done writing initial condition ...'
     ELSE
        CALL output_field("uz",f_buf_i)
@@ -165 +169 @@
          f_u(:), f_W(:), f_uz(:), &  ! IN
          f_uh(:)                         ! OUT
-    REAL(rstd),POINTER :: geopot(:,:), ps(:), rhodz(:,:), u(:,:), W(:,:), uh(:,:), uz(:,:)
+    REAL(rstd),POINTER :: geopot(:,:), ps(:), rhodz(:,:), u(:,:), W(:,:), uh(:,:), uz(:,:), F_el(:,:)
     INTEGER :: ind
     
@@ -177 +181 @@
        W = f_W(ind)
        uh  = f_uh(ind)
+       F_el  = f_buf_Fel(ind)
        IF(caldyn_eta==eta_mass) THEN
           ps=f_ps(ind)
@@ -182 +187 @@
        END IF
        uz = f_uz(ind)
-       CALL compute_prognostic_vel_to_horiz(geopot,rhodz,u,W,uh,uz)
+       !$OMP BARRIER
+       CALL compute_prognostic_vel_to_horiz(geopot,rhodz,u,W, F_el, uh,uz)
+       !$OMP BARRIER
     END DO
   END SUBROUTINE progonostic_vel_to_horiz
   
-  SUBROUTINE compute_prognostic_vel_to_horiz(Phi, rhodz, u, W, uh, uz)
+  SUBROUTINE compute_prognostic_vel_to_horiz(Phi, rhodz, u, W, F_el, uh, uz)
     USE omp_para
     REAL(rstd), INTENT(IN) :: Phi(iim*jjm,llm+1)
@@ -220 +227 @@
        END DO
 
+       ! We need a barrier here because we compute F_el above and do a vertical average below
+       !$OMP BARRIER
+
        DO l=ll_begin, ll_end ! compute on k levels (full levels)
           DO ij=ij_begin_ext, ij_end_ext

codes/icosagcm/trunk/src/initial/etat0.f90

@@ -1 +1 @@
 MODULE etat0_mod
   USE icosa
+  USE omp_para
   IMPLICIT NONE         
   PRIVATE
@@ -10 +11 @@
 
     PUBLIC :: etat0, init_etat0, etat0_type
+
+! Important notes for OpenMP
+! When etat0 is called, vertical OpenMP parallelism is deactivated. 
+! Therefore only the omp_level_master thread must work, i.e. :
+!   !$OMP BARRIER
+!    DO ind=1,ndomain
+!      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
+!      ...
+!    END DO
+!   !$OMP BARRIER
+! There MUST be NO OMP BARRIER inside the DO-LOOP or any routine it calls.
 
 CONTAINS
   
   SUBROUTINE init_etat0
-  USE etat0_database_mod, ONLY: init_etat0_database => init_etat0 
-  USE etat0_start_file_mod, ONLY: init_etat0_start_file => init_etat0 
-  USE etat0_heldsz_mod, ONLY: init_etat0_held_suarez => init_etat0 
-  IMPLICIT NONE
-
+    USE etat0_database_mod, ONLY: init_etat0_database => init_etat0 
+    USE etat0_start_file_mod, ONLY: init_etat0_start_file => init_etat0 
+    USE etat0_heldsz_mod, ONLY: init_etat0_held_suarez => init_etat0 
+    
     CALL getin("etat0",etat0_type)
 
@@ -75 +86 @@
     USE etat0_start_file_mod, ONLY : etat0_start_file=>etat0  
 
-    IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
     TYPE(t_field),POINTER :: f_mass(:)
@@ -160 +170 @@
     END SELECT
 
-!       !$OMP BARRIER
     IF(autoinit_mass) THEN
        DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
+          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
           CALL swap_dimensions(ind)
           CALL swap_geometry(ind)
@@ -175 +184 @@
   SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
     USE theta2theta_rhodz_mod
-    IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
     TYPE(t_field),POINTER :: f_mass(:)
@@ -200 +208 @@
 
     DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
+      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
+!      IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
@@ -232 +241 @@
     USE exner_mod
     USE omp_para
-    IMPLICIT NONE
     REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
     REAL(rstd),INTENT(IN)  :: temp(iim*jjm,llm)
@@ -248 +256 @@
     CALL compute_pression(ps,p,offset)
     ! flush p
-    !$OMP BARRIER
     DO    l    = ll_begin, ll_end
        DO j=jj_begin-offset,jj_end+offset
@@ -278 +285 @@
        ENDDO
     ENDDO
-    !$OMP BARRIER  
   END SUBROUTINE compute_temperature2entropy
 
@@ -296 +302 @@
     USE etat0_dcmip2016_cyclone_mod, ONLY : compute_dcmip2016_cyclone => compute_etat0
     USE etat0_dcmip2016_supercell_mod, ONLY : compute_dcmip2016_supercell => compute_etat0
-    IMPLICIT NONE
     REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
     REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm)
@@ -328 +333 @@
     autoinit_NH = .NOT. hydrostatic
     w(:,:) = 0
-
-    !$OMP BARRIER
 
     SELECT CASE (TRIM(etat0_type))
@@ -391 +394 @@
     END IF
 
-    !$OMP BARRIER
-
     CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u)
 
@@ -405 +406 @@
 
   SUBROUTINE compute_etat0_isothermal(ngrid, phis, ps, temp, ulon, ulat, q)
-    IMPLICIT NONE  
     INTEGER, INTENT(IN)    :: ngrid
     REAL(rstd),INTENT(OUT) :: phis(ngrid)

codes/icosagcm/trunk/src/output/write_etat0.f90

@@ -5 +5 @@
 CONTAINS
 
-  SUBROUTINE write_etat0(it,f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+  SUBROUTINE write_etat0(it,f_ps,f_phis,f_theta_rhodz,f_u, f_q, f_geopot, f_W)
   USE icosa
   USE restart_mod
@@ -20 +20 @@
     TYPE(t_field),POINTER :: f_u(:)
     TYPE(t_field),POINTER :: f_q(:)
+    TYPE(t_field),POINTER, OPTIONAL :: f_geopot(:)
+    TYPE(t_field),POINTER, OPTIONAL :: f_W(:)
   
     TYPE(t_field),POINTER,SAVE :: f_ulon(:)
@@ -41 +43 @@
     CALL transfert_request(f_u,req_e1_vect)
     CALL un2ulonlat(f_u, f_ulon, f_ulat)
-    
-    CALL write_restart(it,f_ps,f_phis,f_theta_rhodz_1d,f_u, f_ulon, f_ulat, f_q)
 
+    IF(hydrostatic) THEN
+       CALL write_restart(it,f_ps,f_phis,f_theta_rhodz_1d,f_u, f_ulon, f_ulat, f_q)
+    ELSE
+       CALL write_restart(it,f_ps,f_phis,f_theta_rhodz_1d,f_u, f_ulon, f_ulat, f_q, f_geopot, f_W)
+    END IF
     CALL deallocate_field(f_ulon)
     CALL deallocate_field(f_ulat)

codes/icosagcm/trunk/src/time/timeloop_gcm.f90

@@ -219 +219 @@
     END IF
     
-    IF (write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q)  ! FIXME : write_start undefined
+!    IF (write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q)
+    IF (write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q, f_geopot, f_W)
    
     CALL write_output_fields_basic(.TRUE., f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
@@ -337 +338 @@
     END DO
     
-    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q) 
+!    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q) 
+    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q, f_geopot, f_W) 
     
     CALL check_conserve_detailed(it, AAM_dyn, &