Changeset 326 for codes


Ignore:
Timestamp:
02/02/15 01:12:53 (9 years ago)
Author:
dubos
Message:

Detailed energy/AAM diagnostics

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

Legend:

Unmodified
Added
Removed

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

    r295 r326  
    77    TYPE(t_field),POINTER,SAVE :: f_vort(:) 
    88    TYPE(t_field),POINTER,SAVE :: f_rhodz(:)  
    9     
    10   PUBLIC init_check_conserve, check_conserve  
     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 
    1113    REAL(rstd),SAVE:: mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0 
    12 !$OMP THREADPRIVATE(mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0) 
     14!$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0) 
    1315      
    1416CONTAINS  
    … …  
    1719 
    1820  SUBROUTINE init_check_conserve 
    19   USE icosa 
    20   IMPLICIT NONE 
     21    USE icosa 
     22    USE getin_mod 
     23    USE omp_para, ONLY : is_master 
     24    IMPLICIT NONE 
     25    CHARACTER(LEN=255) :: check_type_str 
    2126    CALL allocate_field(f_pk,field_t,type_real,llm) 
    2227    CALL allocate_field(f_p,field_t,type_real,llm+1) 
    … …  
    2429    CALL allocate_field(f_rhodz,field_t,type_real,llm) 
    2530    CALL allocate_field(f_vort,field_z,type_real,llm) 
     31     
     32    check_type_str='basic' 
     33    CALL getin("check_conservation",check_type_str) 
     34    SELECT CASE(TRIM(check_type_str)) 
     35    CASE('basic') 
     36       check_type=check_basic 
     37    CASE('detailed') 
     38       check_type=check_detailed 
     39    CASE DEFAULT 
     40       IF(is_master) PRINT*,'Bad selector <', TRIM(check_type_str), & 
     41            '> for variable check_conservation : options are <basic>, <detailed>' 
     42       STOP 
     43    END SELECT 
    2644  END SUBROUTINE init_check_conserve 
     45 
     46!--------------------------------- Basic check -------------------------------- 
    2747 
    2848  SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it)  
    … …  
    4060    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    4161    TYPE(t_field),POINTER :: f_phis(:) 
    42     INTEGER::it 
     62    INTEGER, INTENT(IN) :: it 
    4363 
    4464    REAL(rstd),POINTER :: p(:,:),rhodz(:,:)  
    4565    INTEGER::ind,ierr 
    46     REAL(rstd) :: mtot, rmsdpdt 
    47     REAL(rstd) :: etot, stot, angtot, rmsvtot, ztot 
     66    REAL(rstd) :: mtot, angtot, rmsdpdt 
     67    REAL(rstd) :: etot, stot, ang_mass, ang_vit, rmsvtot, ztot 
    4868 
    4969    CALL transfert_request(f_ue,req_e1_vect) 
    … …  
    6383    CALL check_PV(ztot)  
    6484    CALL exner(f_ps,f_p,f_pks,f_pk) 
    65     CALL check_EN(f_ue,f_theta_rhodz,f_phis, etot, stot, angtot, rmsvtot)  
     85    CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vit, rmsvtot)  
     86    angtot  = ang_mass + ang_vit  
    6687 
    6788    IF (is_master) THEN  
    … …  
    7091          mtot0 = mtot 
    7192          etot0 = etot  
    72           angtot0  = angtot  
     93          angtot0  = angtot 
    7394          stot0 = stot  
    7495       END IF 
    … …  
    92113    END IF 
    93114  END SUBROUTINE check_conserve 
    94    
     115 
     116!------------------------ Detailed check (Lauritzen et al., 2014) ----------------------------- 
     117 
     118  SUBROUTINE check_conserve_detailed(tag, f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it)  
     119    USE icosa  
     120    USE pression_mod 
     121    USE vorticity_mod 
     122    USE caldyn_gcm_mod 
     123    USE exner_mod  
     124    USE mpipara, ONLY : is_mpi_root, comm_icosa 
     125    IMPLICIT NONE 
     126    CHARACTER(LEN=*), INTENT(IN) :: tag 
     127    TYPE(t_field),POINTER :: f_ps(:) 
     128    TYPE(t_field),POINTER :: f_dps(:) 
     129    TYPE(t_field),POINTER :: f_ue(:) 
     130    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
     131    TYPE(t_field),POINTER :: f_phis(:) 
     132    INTEGER, INTENT(IN) ::it 
     133 
     134    REAL(rstd),POINTER :: p(:,:),rhodz(:,:)  
     135    INTEGER::ind,ierr 
     136    REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, AAM_mass, AAM_vel 
     137     
     138    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  
     142 
     143       CALL transfert_request(f_ue,req_e1_vect) 
     144       CALL pression(f_ps,f_p) 
     145       DO ind=1,ndomain  
     146          IF (.NOT. assigned_domain(ind)) CYCLE 
     147          CALL swap_dimensions(ind) 
     148          CALL swap_geometry(ind) 
     149          p=f_p(ind)  
     150          rhodz=f_rhodz(ind) 
     151          CALL compute_rhodz(p,rhodz)  
     152       END DO 
     153       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)  
     155        
     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  
     168       END IF 
     169    END IF 
     170  END SUBROUTINE check_conserve_detailed 
     171  
    95172!--------------------------------------------------------------------- 
     173 
     174  SUBROUTINE global_sum(loc_sum, glob_sum) 
     175    USE mpi_mod 
     176    USE mpipara 
     177    USE transfert_omp_mod 
     178    IMPLICIT NONE 
     179    REAL(rstd), INTENT(IN) :: loc_sum 
     180    REAL(rstd), INTENT(OUT) :: glob_sum 
     181    REAL(rstd) :: sum 
     182 
     183    CALL reduce_sum_omp(loc_sum, sum) 
     184    IF (using_mpi) THEN 
     185!$OMP MASTER     
     186      CALL MPI_REDUCE(sum, glob_sum, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 
     187!$OMP END MASTER 
     188    ELSE 
     189      glob_sum=loc_sum 
     190    ENDIF  
     191  END SUBROUTINE global_sum 
    96192 
    97193  SUBROUTINE check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt) 
    … …  
    128224      ENDDO 
    129225    ENDDO 
    130  
    131     IF (using_mpi) THEN 
    132         CALL reduce_sum_omp(mloc, mloc_mpi) 
    133         CALL reduce_sum_omp(rmsloc, rmsloc_mpi) 
    134 !$OMP MASTER     
    135       CALL MPI_REDUCE(mloc_mpi, mtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 
    136       CALL MPI_REDUCE(rmsloc_mpi, rmsdpdt, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 
    137 !$OMP END MASTER 
    138     ELSE 
    139       mtot=mloc 
    140       rmsdpdt=rmsloc 
    141     ENDIF 
    142      
     226     
     227    CALL global_sum(mloc, mtot) 
     228    CALL global_sum(rmsloc, rmsdpdt) 
    143229  END SUBROUTINE check_mass_conserve 
    144230 
    145231!--------------------------------------------------------------------- 
    146232 
    147   SUBROUTINE check_en(f_ue,f_theta_rhodz,f_phis, etot, stot, angtot, rmsvtot)   
     233  SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, & 
     234                      stot, AAM_mass_tot, AAM_vel_tot, rmsvtot)   
    148235  USE icosa 
    149236  USE pression_mod 
    150237  USE vorticity_mod 
    151   USE mpi_mod 
    152   USE mpipara 
    153   USE transfert_omp_mod 
    154238  IMPLICIT NONE 
    155239    TYPE(t_field), POINTER :: f_ue(:) 
    156240    TYPE(t_field), POINTER :: f_theta_rhodz(:) 
    157241    TYPE(t_field), POINTER :: f_phis(:) 
    158     REAL(rstd),INTENT(OUT) :: etot, stot, angtot, rmsvtot 
     242    REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, rmsvtot 
    159243    REAL(rstd), POINTER :: ue(:,:) 
    160244    REAL(rstd), POINTER :: p(:,:)  
    … …  
    163247    REAL(rstd), POINTER :: phis(:)  
    164248    REAL(rstd), POINTER :: rhodz(:,:)  
    165     REAL(rstd) :: e, s, ang, rmsv 
    166     REAL(rstd) :: e_mpi, s_mpi, ang_mpi, rmsv_mpi 
     249    REAL(rstd) :: e, s, AAM_mass, AAM_vel, rmsv 
    167250    INTEGER :: ind 
    168251 
    169     e = 0 ; s = 0 ; ang = 0 ; rmsv = 0 
     252    e = 0 ; s = 0 ; AAM_mass = 0 ; AAM_vel=0 ; rmsv = 0 
    170253 
    171254    DO ind=1,ndomain  
    … …  
    179262      pk=f_pk(ind)  
    180263      phis=f_phis(ind)  
    181       CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis, e, s, ang, rmsv) 
     264      CALL compute_energy(ind,ue,p,rhodz,theta_rhodz,pk,phis, & 
     265           e, s, AAM_mass, AAM_vel, rmsv) 
    182266    END DO  
    183267 
    184     IF (using_mpi) THEN 
    185       CALL reduce_sum_omp(e, e_mpi) 
    186       CALL reduce_sum_omp(s, s_mpi) 
    187       CALL reduce_sum_omp(ang, ang_mpi) 
    188       CALL reduce_sum_omp(rmsv, rmsv_mpi) 
    189 !$OMP MASTER     
    190       CALL MPI_REDUCE(e_mpi, etot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 
    191       CALL MPI_REDUCE(s_mpi, stot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 
    192       CALL MPI_REDUCE(ang_mpi, angtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 
    193       CALL MPI_REDUCE(rmsv_mpi, rmsvtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 
    194 !$OMP END MASTER 
    195     ELSE 
    196       etot=e 
    197       stot=s 
    198       angtot=ang 
    199       rmsvtot=rmsv 
    200     ENDIF  
    201  
    202   END SUBROUTINE check_en  
    203  
    204   SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, ang, rmsv) 
    205   USE icosa 
    206   USE disvert_mod 
    207   USE wind_mod 
    208   USE omp_para 
    209   IMPLICIT NONE 
     268    CALL global_sum(e, etot) 
     269    CALL global_sum(s, stot) 
     270    CALL global_sum(rmsv, rmsvtot) 
     271    CALL global_sum(AAM_mass, AAM_mass_tot) 
     272    CALL global_sum(AAM_vel, AAM_vel_tot) 
     273  END SUBROUTINE check_energy 
     274 
     275  SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, rmsv) 
     276    USE icosa 
     277    USE disvert_mod 
     278    USE wind_mod 
     279    USE omp_para 
     280    IMPLICIT NONE 
    210281    INTEGER,INTENT(IN)::ind  
    211282    REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm) 
    … …  
    215286    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm) 
    216287    REAL(rstd),INTENT(IN) :: phis(iim*jjm) 
    217     REAL(rstd),INTENT(INOUT) :: e 
    218     REAL(rstd),INTENT(INOUT) :: s 
    219     REAL(rstd),INTENT(INOUT) :: ang 
    220     REAL(rstd),INTENT(INOUT) :: rmsv 
    221      
    222     REAL(rstd):: theta(iim*jjm,llm) 
    223     REAL(rstd)::KE(iim*jjm,llm)  
     288    REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, rmsv 
     289     
    224290    REAL(rstd) :: ucenter(iim*jjm,3,llm) 
    225291    REAL(rstd) :: ulon(iim*jjm,llm) 
    226292    REAL(rstd) :: ulat(iim*jjm,llm) 
    227  
    228     REAL(rstd)::masse(iim*jjm,llm)   
    229     REAL(rstd)::rad,radomeg,lat,lon 
    230     REAL(rstd) ::xyz(3) 
    231     INTEGER :: i,j,ij,l,ij2 
    232  
    233    
    234     DO l = ll_begin, ll_end 
    235       DO j=jj_begin,jj_end 
    236         DO i=ii_begin,ii_end 
    237           ij=(j-1)*iim+i 
    238           masse(ij,l) = rhodz(ij,l)*Ai(ij)  
    239           theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    240           IF (domain(ind)%own(i,j)) THEN 
    241             s = s + Ai(ij)*theta_rhodz(ij,l) 
    242           END IF  
    243         END DO  
    244       END DO 
    245     END DO  
    246  
    247     DO l = ll_begin,ll_end 
    248       DO j=jj_begin,jj_end 
    249         DO i=ii_begin,ii_end 
    250           ij=(j-1)*iim+i 
    251           IF (domain(ind)%own(i,j)) THEN 
    252             KE(ij,l)=   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   & 
    253                                       le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           & 
    254                                       le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           & 
    255                                       le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        & 
    256                                       le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     & 
    257                                       le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    258             rmsv = rmsv + 2*masse(ij,l)*KE(ij,l)  
    259           END IF  
    260         END DO 
    261       END DO  
    262     END DO   
    263  
    264     DO l = ll_begin,ll_end 
    265       DO j=jj_begin,jj_end 
    266         DO i=ii_begin,ii_end 
    267           ij=(j-1)*iim+i 
    268           IF (domain(ind)%own(i,j)) THEN 
    269             e = e + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l)) 
    270           END IF  
    271         END DO  
    272       END DO  
    273     END DO  
    274  
    275 !------------------------------ ANGULAR VELOCITY  
     293     
     294    REAL(rstd) :: mass_ij, theta_ij, KE_ij, rad,radomeg,lat 
     295    INTEGER :: i,j,ij,l 
     296     
    276297    CALL compute_wind_centered(u,ucenter) 
    277298    CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat) 
    278299    rad = radius  
    279300    radomeg = rad*omega 
    280  
     301     
    281302    DO l = ll_begin,ll_end 
    282       DO j=jj_begin,jj_end 
    283         DO i=ii_begin,ii_end 
    284           ij=(j-1)*iim+i 
    285           IF (domain(ind)%own(i,j)) THEN 
    286             lat=lat_i(ij)   
    287             ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat)) 
    288           END IF  
    289         END DO 
    290       END DO 
    291     END DO     
    292            
    293   END SUBROUTINE compute_en 
     303       DO j=jj_begin,jj_end 
     304          DO i=ii_begin,ii_end 
     305             ij=(j-1)*iim+i 
     306             IF (domain(ind)%own(i,j)) THEN 
     307                lat = lat_i(ij)   
     308                mass_ij = rhodz(ij,l)*Ai(ij)  
     309                theta_ij = theta_rhodz(ij,l)/rhodz(ij,l) 
     310                KE_ij = 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   & 
     311                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           & 
     312                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           & 
     313                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        & 
     314                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     & 
     315                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
     316                rmsv = rmsv + 2*mass_ij*KE_ij 
     317                s    = s + Ai(ij)*theta_rhodz(ij,l) 
     318                e    = e + mass_ij*(phis(ij)+theta_ij*pk(ij,l)+KE_ij) 
     319                AAM_mass = AAM_mass + rad*cos(lat)*mass_ij*radomeg*cos(lat) 
     320                AAM_vel  = AAM_vel  + rad*cos(lat)*mass_ij*ulon(ij,l) 
     321             END IF 
     322          END DO 
     323       END DO 
     324    END DO 
     325     
     326END SUBROUTINE compute_energy 
    294327 
    295328!--------------------------------------------------------------------- 

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

    r295 r326  
    211211   
    212212  DO it=itau0+1,itau0+itaumax 
    213      
     213      
     214    CALL check_conserve_detailed('detailed_budget 0', & 
     215         f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
     216 
    214217    IF (xios_output) CALL xios_update_calendar(it) 
    215218    IF (it==itau0+1 .OR. MOD(it,itau_sync)==0) THEN 
    … …  
    256259       END SELECT 
    257260    END DO 
     261 
     262    CALL check_conserve_detailed('detailed_budget 1', & 
     263         f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
    258264 
    259265    IF (MOD(it,itau_dissip)==0) THEN 
    … …  
    281287    END IF 
    282288 
     289    CALL check_conserve_detailed('detailed_budget 2', & 
     290         f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
     291 
    283292    IF(MOD(it,itau_adv)==0) THEN 
    284293 
Note: See TracChangeset for help on using the changeset viewer.