Changeset 365 for codes


Ignore:
Timestamp:
10/12/15 11:19:11 (9 years ago)
Author:
dubos
Message:

Cleanup detailed AAM budget

Location:
codes/icosagcm/trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r326 r365  
    33  IMPLICIT NONE  
    44 
    5   PRIVATE  
    6     TYPE(t_field),POINTER,SAVE :: f_pk(:), f_pks(:), f_p(:) 
    7     TYPE(t_field),POINTER,SAVE :: f_vort(:) 
    8     TYPE(t_field),POINTER,SAVE :: f_rhodz(:)  
    9  
    10     INTEGER, PARAMETER :: check_basic=1, check_detailed=2 
    11     INTEGER, SAVE :: check_type 
    12   PUBLIC init_check_conserve, check_conserve, check_conserve_detailed 
    13     REAL(rstd),SAVE:: mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0 
    14 !$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0) 
    15       
     5  PRIVATE 
     6   
     7  INTEGER, PARAMETER, PUBLIC :: check_basic=1, check_detailed=2, & 
     8       AAM_dissip=1, AAM_dyn=2, AAM_phys=3 
     9  INTEGER, SAVE :: check_type 
     10 
     11  TYPE(t_field),POINTER,SAVE :: f_pk(:), f_pks(:), f_p(:) 
     12  TYPE(t_field),POINTER,SAVE :: f_vort(:) 
     13  TYPE(t_field),POINTER,SAVE :: f_rhodz(:)  
     14   
     15  REAL(rstd),SAVE :: AAM_mass, AAM_mass_old, AAM_vel, AAM_vel_old, & 
     16       AAM_mass_source(3), AAM_vel_source(3) ! read/written only IF is_master 
     17  REAL(rstd),SAVE :: mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0 
     18  !$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0) 
     19   
     20  PUBLIC :: init_check_conserve, check_conserve_detailed, check_conserve 
    1621CONTAINS  
    1722 
    … …  
    1924 
    2025  SUBROUTINE init_check_conserve 
    21     USE icosa 
    2226    USE getin_mod 
    2327    USE omp_para, ONLY : is_master 
    24     IMPLICIT NONE 
    2528    CHARACTER(LEN=255) :: check_type_str 
    2629    CALL allocate_field(f_pk,field_t,type_real,llm) 
    … …  
    4750 
    4851  SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it)  
    49   USE icosa  
    5052  USE pression_mod 
    5153  USE vorticity_mod 
    5254  USE caldyn_gcm_mod 
    5355  USE exner_mod  
    54   USE mpipara, ONLY : is_mpi_master, comm_icosa 
    55   USE omp_para 
    56   IMPLICIT NONE 
     56!  USE mpipara, ONLY : is_mpi_master, comm_icosa 
     57  USE omp_para, ONLY : is_master 
    5758    TYPE(t_field),POINTER :: f_ps(:) 
    5859    TYPE(t_field),POINTER :: f_dps(:) 
    … …  
    6364 
    6465    REAL(rstd),POINTER :: p(:,:),rhodz(:,:)  
    65     INTEGER::ind,ierr 
     66    INTEGER :: ind,ierr 
    6667    REAL(rstd) :: mtot, angtot, rmsdpdt 
    67     REAL(rstd) :: etot, stot, ang_mass, ang_vit, rmsvtot, ztot 
     68    REAL(rstd) :: etot, stot, ang_mass, ang_vel, rmsvtot, ztot 
    6869 
    6970    CALL transfert_request(f_ue,req_e1_vect) 
    … …  
    8384    CALL check_PV(ztot)  
    8485    CALL exner(f_ps,f_p,f_pks,f_pk) 
    85     CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vit, rmsvtot)  
    86     angtot  = ang_mass + ang_vit  
    87  
    88     IF (is_master) THEN  
     86    CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, rmsvtot)  
     87 
     88    IF (is_master) THEN ! is_master = is_omp_master && is_mpi_master 
     89       ! AAM_mass and AAM_vel must be read/written only IF is_master 
     90       AAM_mass = ang_mass 
     91       AAM_vel = ang_vel 
     92       angtot  = ang_mass + ang_vel 
    8993       IF ( it == itau0  ) THEN  
    9094          ztot0 = ztot 
    … …  
    9397          angtot0  = angtot 
    9498          stot0 = stot  
     99          AAM_mass_old=AAM_mass 
     100          AAM_vel_old=AAM_vel 
     101          AAM_mass_source=0. 
     102          AAM_vel_source=0. 
    95103       END IF 
    96  
    97104       rmsvtot=SQRT(rmsvtot/mtot)  
    98105       ztot=ztot/ztot0-1. ; mtot=mtot/mtot0-1.  
    … …  
    107114        
    1081154000   FORMAT(10x,'masse',5x,'rmsdpdt',5x,'energie',5x,'enstrophie'  & 
    109             ,5x,'entropie',5x,'rmsv',5x,'mt.ang',/,'GLOB  '                & 
     116            ,5x,'entropie',5x,'rmsv',5x,'mt.ang',/,'GLOB  '          & 
    110117            ,e10.3,e13.6,5e13.3/)      
    111        close(134) 
    112        
     118       CLOSE(134) 
     119 
     120       IF(check_type == check_detailed) THEN 
     121          !ang_mass = AAM_mass-AAM_mass_old  
     122          !WRITE(*,'(A,3E12.4)') 'AAM_mass sanity check', SUM(AAM_mass_source), ang_mass, SUM(AAM_mass_source)-ang_mass 
     123          !ang_vel = AAM_vel-AAM_vel_old  
     124          !WRITE(*,'(A,3E12.4)') 'AAM_vel  sanity check', SUM(AAM_vel_source), ang_vel, SUM(AAM_vel_source)-ang_vel 
     125          WRITE(*,'(A,6E12.4)') 'AAM_mass : time,old,new,dissip,dyn,phys', dt*it, AAM_mass_old, AAM_mass, & 
     126               AAM_mass_source(AAM_dissip), AAM_mass_source(AAM_dyn), AAM_mass_source(AAM_phys) 
     127          WRITE(*,'(A,6E12.4)') 'AAM_vel  : time,old,new,dissip,dyn,phys', dt*it, AAM_vel_old, AAM_vel, & 
     128               AAM_vel_source(AAM_dissip), AAM_vel_source(AAM_dyn), AAM_vel_source(AAM_phys) 
     129          AAM_mass_old=AAM_mass 
     130          AAM_vel_old=AAM_vel 
     131          AAM_mass_source=0. 
     132          AAM_vel_source=0. 
     133       END IF 
    113134    END IF 
     135 
    114136  END SUBROUTINE check_conserve 
    115137 
    116138!------------------------ Detailed check (Lauritzen et al., 2014) ----------------------------- 
    117139 
    118   SUBROUTINE check_conserve_detailed(tag, f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it)  
    119     USE icosa  
     140  SUBROUTINE check_conserve_detailed(it,tag, f_ps,f_dps,f_ue,f_theta_rhodz,f_phis)  
    120141    USE pression_mod 
    121142    USE vorticity_mod 
    122143    USE caldyn_gcm_mod 
    123144    USE exner_mod  
    124     USE mpipara, ONLY : is_mpi_root, comm_icosa 
    125     IMPLICIT NONE 
    126     CHARACTER(LEN=*), INTENT(IN) :: tag 
     145!    USE mpipara, ONLY : is_mpi_root, comm_icosa 
     146    USE omp_para, ONLY : is_master 
     147    INTEGER, INTENT(IN) :: tag 
    127148    TYPE(t_field),POINTER :: f_ps(:) 
    128149    TYPE(t_field),POINTER :: f_dps(:) 
    … …  
    134155    REAL(rstd),POINTER :: p(:,:),rhodz(:,:)  
    135156    INTEGER::ind,ierr 
    136     REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, AAM_mass, AAM_vel 
     157    REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, ang_mass, ang_vel 
    137158     
    138159    IF(check_type == check_detailed) THEN 
    139        etot=0.0; stot=0.0; rmsv=0.0  
    140        AAM_vel=0. ; AAM_vel=0. 
    141        ztot = 0.0  
    142160 
    143161       CALL transfert_request(f_ue,req_e1_vect) 
    … …  
    152170       END DO 
    153171       CALL exner(f_ps,f_p,f_pks,f_pk) 
    154        CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, AAM_mass, AAM_vel, rmsv)  
     172       CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, rmsv)  
    155173        
    156        IF (is_mpi_root) THEN  
    157           !$OMP MASTER        
    158           IF ( it == itau0  ) Then  
    159              ztot0 = ztot 
    160              mtot0 = mtot 
    161              etot0 = etot  
    162              angtot0  = AAM_mass + AAM_vel 
    163              stot0 = stot  
    164           END IF 
    165            
    166           PRINT *, TRIM(tag), AAM_mass, AAM_vel 
    167           !$OMP END MASTER  
     174       IF (is_master) THEN 
     175          ! XXX=mass,vel 
     176          ! AAM_XXX = value of total AAM_XXX at last check 
     177          ! ang_XXX = value of total AAM_XXX now 
     178          ! AAM_XXX_source(tag) = source of AAA_XXX due to tag=dissip,dyn,phys 
     179          !WRITE(*,'(A,3E12.2)') 'AAM_mass', tag, ang_mass, AAM_mass, ang_mass-AAM_mass ! FIXME 
     180          !WRITE(*,'(A,3E12.2)') 'AAM_vel ', tag, ang_vel, AAM_vel, ang_vel-AAM_vel ! FIXME 
     181          AAM_mass_source(tag) = AAM_mass_source(tag) + ang_mass - AAM_mass 
     182          AAM_vel_source(tag) = AAM_vel_source(tag) + ang_vel - AAM_vel 
     183          AAM_mass = ang_mass 
     184          AAM_vel = ang_vel 
    168185       END IF 
    169186    END IF 
    … …  
    176193    USE mpipara 
    177194    USE transfert_omp_mod 
    178     IMPLICIT NONE 
    179195    REAL(rstd), INTENT(IN) :: loc_sum 
    180196    REAL(rstd), INTENT(OUT) :: glob_sum 
    … …  
    188204    ELSE 
    189205      glob_sum=loc_sum 
    190     ENDIF  
     206    ENDIF 
    191207  END SUBROUTINE global_sum 
    192208 
    … …  
    195211  USE mpipara 
    196212  USE transfert_omp_mod 
    197   USE icosa 
    198213  USE omp_para 
    199214  IMPLICIT NONE 
    … …  
    233248  SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, & 
    234249                      stot, AAM_mass_tot, AAM_vel_tot, rmsvtot)   
    235   USE icosa 
    236250  USE pression_mod 
    237251  USE vorticity_mod 
    238   IMPLICIT NONE 
    239252    TYPE(t_field), POINTER :: f_ue(:) 
    240253    TYPE(t_field), POINTER :: f_theta_rhodz(:) 
    … …  
    274287 
    275288  SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, rmsv) 
    276     USE icosa 
    277289    USE disvert_mod 
    278290    USE wind_mod 
    279291    USE omp_para 
    280     IMPLICIT NONE 
    281292    INTEGER,INTENT(IN)::ind  
    282293    REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm) 
    … …  
    324335    END DO 
    325336     
    326 END SUBROUTINE compute_energy 
     337  END SUBROUTINE compute_energy 
    327338 
    328339!--------------------------------------------------------------------- 
    329  
     340   
    330341  SUBROUTINE check_PV(ztot) 
    331   USE icosa 
    332   USE mpi_mod 
    333   USE mpipara 
    334   USE transfert_omp_mod 
    335   IMPLICIT NONE 
    336      REAL(rstd),INTENT(OUT) :: ztot 
    337      REAL(rstd), POINTER :: vort(:,:) 
    338      REAL(rstd), POINTER :: rhodz(:,:)  
    339      INTEGER :: ind 
    340      REAL(rstd) :: z,z_mpi 
    341  
    342      z=0 
    343      DO ind=1,ndomain  
     342    REAL(rstd),INTENT(OUT) :: ztot 
     343    REAL(rstd), POINTER :: vort(:,:) 
     344    REAL(rstd), POINTER :: rhodz(:,:)  
     345    INTEGER :: ind 
     346    REAL(rstd) :: z,z_mpi 
     347     
     348    z=0 
     349    DO ind=1,ndomain  
    344350       IF (.NOT. assigned_domain(ind)) CYCLE 
    345351       CALL swap_dimensions(ind) 
    … …  
    348354       rhodz=f_rhodz(ind)  
    349355       CALL compute_PV(vort,rhodz,z)  
    350      ENDDO 
    351  
    352     IF (using_mpi) THEN 
    353       CALL reduce_sum_omp(z, z_mpi) 
    354 !$OMP MASTER     
    355       CALL MPI_REDUCE(z_mpi, ztot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 
    356 !$OMP END MASTER 
    357     ELSE 
    358       ztot=z 
    359     ENDIF  
    360   
     356    ENDDO 
     357    CALL global_sum(z, ztot) 
    361358  END SUBROUTINE check_PV  
    362359 
    363360  SUBROUTINE compute_PV(vort,rhodz,z)  
    364   USE icosa 
    365361  USE disvert_mod 
    366362  USE omp_para 
    367   IMPLICIT NONE  
    368363    REAL(rstd),INTENT(IN) :: vort(2*iim*jjm,llm) 
    369364    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)  
    … …  
    402397 
    403398   SUBROUTINE compute_rhodz(p,rhodz)  
    404    USE icosa  
    405    IMPLICIT NONE  
    406399     REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1) 
    407400     REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm) 

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

    r364 r365  
    225225 
    226226       IF (mod(it,itau_out)==0 ) THEN 
     227          CALL transfert_request(f_u,req_e1_vect) 
    227228          CALL write_output_fields_basic(f_ps, f_u, f_q) 
    228229       ENDIF 
    229  
    230        CALL check_conserve_detailed('detailed_budget 0', & 
    231             f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
    232230 
    233231       CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 
    … …  
    239237          CALL HEVI_scheme(it, fluxt_zero) 
    240238       END SELECT 
    241         
    242        CALL check_conserve_detailed('detailed_budget 1', & 
    243             f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
    244239        
    245240       IF (MOD(it,itau_dissip)==0) THEN 
    … …  
    257252          ENDIF 
    258253           
     254          CALL check_conserve_detailed(it, AAM_dyn, & 
     255               f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 
     256        
    259257          CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz) 
    260258           
    … …  
    264262             CALL euler_scheme(.FALSE.)  ! update only u, theta 
    265263          END IF 
    266        END IF 
    267         
    268        CALL check_conserve_detailed('detailed_budget 2', & 
    269             f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
     264 
     265          CALL check_conserve_detailed(it, AAM_dissip, & 
     266               f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 
     267       END IF 
    270268        
    271269       IF(MOD(it,itau_adv)==0) THEN 
    … …  
    284282       END IF 
    285283        
    286        IF (MOD(it,itau_check_conserv)==0) THEN 
    287           CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
    288        ENDIF 
    289         
    290284       IF (MOD(it,itau_physics)==0) THEN 
     285          CALL check_conserve_detailed(it, AAM_dyn, & 
     286            f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 
    291287          CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)         
     288          CALL check_conserve_detailed(it, AAM_phys, & 
     289               f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 
    292290          !$OMP MASTER 
    293291          IF (first_physic) CALL SYSTEM_CLOCK(start_clock) 
    … …  
    295293          first_physic=.FALSE. 
    296294       END IF 
    297         
     295 
     296       IF (MOD(it,itau_check_conserv)==0) THEN 
     297          CALL check_conserve_detailed(it, AAM_dyn, & 
     298               f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 
     299          CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
     300       ENDIF        
    298301    END DO 
    299302     
    300303    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q)  
    301304     
     305    CALL check_conserve_detailed(it, AAM_dyn, & 
     306         f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 
    302307    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)   
    303308     
    … …  
    316321    INTEGER :: it, itau0, itaumax, start_clock, stop_clock, rate_clock, throughput 
    317322    REAL :: per_step,total, elapsed 
    318     WRITE(*,'(A,I7,A,F8.1)') "It No :",it,"   t :",dt*it 
     323    WRITE(*,'(A,I7,A,F14.1)') "It No :",it,"   t :",dt*it 
    319324    IF(MOD(it,10)==0) THEN 
    320325       CALL SYSTEM_CLOCK(stop_clock) 
Note: See TracChangeset for help on using the changeset viewer.