Changeset 428 for codes


Ignore:
Timestamp:
06/14/16 21:54:26 (8 years ago)
Author:
dubos
Message:

theta-related cleanup

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

Legend:

Unmodified
Added
Removed

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

    r362 r428  
    3333    USE icosa 
    3434    USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass 
    35     USE exner_mod 
    3635    USE trace 
    3736    USE omp_para 
    … …  
    121120    USE icosa 
    122121    USE disvert_mod 
    123     USE exner_mod 
    124122    USE trace 
    125123    USE omp_para 

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

    r295 r428  
    11MODULE etat0_academic_mod 
    2  
    3  
    4  
    5  
     2  USE icosa 
     3  IMPLICIT NONE 
     4 
     5  PRIVATE 
     6   
     7  PUBLIC :: etat0 
     8   
    69CONTAINS 
    710   
    … …  
    912  USE icosa 
    1013  USE kinetic_mod 
    11   IMPLICIT NONE 
    1214    TYPE(t_field),POINTER,SAVE :: f_ps(:) 
    1315    TYPE(t_field),POINTER,SAVE :: f_phis(:) 
    … …  
    2123    REAL(rstd),POINTER :: temp(:) 
    2224    INTEGER :: ind 
    23      
    24      
     25         
    2526    CALL allocate_field(f_ps,field_t,type_real) 
    2627    CALL allocate_field(f_phis,field_t,type_real) 
    … …  
    4041     
    4142  END SUBROUTINE test_etat0_academic 
    42     
    4343     
    4444  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    4545  USE icosa 
    46   IMPLICIT NONE 
    4746    TYPE(t_field),POINTER :: f_ps(:) 
    4847    TYPE(t_field),POINTER :: f_phis(:) 
    … …  
    5655    REAL(rstd),POINTER :: u(:,:) 
    5756    INTEGER :: ind 
    58      
     57 
     58    PRINT *, 'etat0_academic needs an upgrade for 4D theta' 
     59    PRINT *, 'STOP in etat0_academic.f90/etat0' 
     60    STOP 
     61 
    5962    DO ind=1,ndomain 
    6063      IF (.NOT. assigned_domain(ind)) CYCLE 
    … …  
    7780  USE geopotential_mod 
    7881  USE theta2theta_rhodz_mod 
    79   IMPLICIT NONE   
    8082  REAL(rstd),INTENT(OUT) :: ps(iim*jjm) 
    8183  REAL(rstd),INTENT(OUT) :: phis(iim*jjm) 
    … …  
    99101  REAL(rstd) :: fact(3*iim*jjm) 
    100102  REAL(rstd) :: ut(3*iim*jjm,llm) 
    101    
    102    
     103     
    103104    DO l=1,llm 
    104105       zsig=ap(l)/preff+bp(l) 
    … …  
    137138    CALL compute_pression(ps,p,1)      
    138139!$OMP BARRIER 
    139     CALL compute_exner(ps,p,pks,pk,1)   
    140 !$OMP BARRIER 
    141     CALL compute_geopotential(phis,pks,pk,theta,phi,1) 
     140    CALL compute_geopotential(phis,ps,theta,phi,1) 
    142141 
    143142 

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

    r359 r428  
    2323 
    2424  SUBROUTINE test_etat0_heldsz 
    25     USE icosa 
    2625    USE kinetic_mod 
    27     IMPLICIT NONE 
    2826    TYPE(t_field),POINTER :: f_ps(:) 
    2927    TYPE(t_field),POINTER :: f_phis(:) 
    … …  
    5048 
    5149  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    52     USE icosa 
    5350    USE theta2theta_rhodz_mod 
    5451    IMPLICIT NONE 
    … …  
    9996 
    10097  SUBROUTINE init_Teq 
    101     USE icosa 
    10298    USE disvert_mod, ONLY : ap,bp 
    103     IMPLICIT NONE 
    10499    REAL(rstd),POINTER :: clat(:)  
    105100    REAL(rstd),POINTER :: theta_eq(:,:)  
    … …  
    159154 
    160155  SUBROUTINE compute_Teq(lat,theta_eq) 
    161     USE icosa 
    162156    USE disvert_mod 
    163     IMPLICIT NONE   
    164157    REAL(rstd),INTENT(IN) :: lat(iim*jjm) 
    165158    REAL(rstd),INTENT(OUT) :: theta_eq(iim*jjm,llm)  
    … …  
    184177 
    185178  SUBROUTINE compute_etat0_heldsz(theta_eq, theta) 
    186     USE icosa 
    187179    USE disvert_mod 
    188     USE pression_mod 
    189     USE exner_mod 
    190     USE geopotential_mod 
    191     USE theta2theta_rhodz_mod 
    192     IMPLICIT NONE   
    193180    REAL(rstd),INTENT(IN) :: theta_eq(iim*jjm,llm)  
    194181    REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)  
    … …  
    211198 
    212199  SUBROUTINE held_suarez(f_ps,f_theta_rhodz,f_u)  
    213     USE icosa 
    214     IMPLICIT NONE 
    215200    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    216201    TYPE(t_field),POINTER :: f_u(:) 
    … …  
    238223 
    239224  SUBROUTINE compute_heldsz(ps,theta_eq,lat, theta_rhodz,u, theta)  
    240     USE icosa 
    241225    USE theta2theta_rhodz_mod 
    242     IMPLICIT NONE 
    243226    REAL(rstd),INTENT(IN)    :: ps(iim*jjm)  
    244227    REAL(rstd),INTENT(IN)    :: theta_eq(iim*jjm,llm)  

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

    r295 r428  
    1111  USE icosa 
    1212  IMPLICIT NONE 
    13     TYPE(t_field), POINTER :: f_ps(:) 
    14     TYPE(t_field), POINTER :: f_p(:) 
    15     TYPE(t_field), POINTER :: f_pks(:) 
    16     TYPE(t_field), POINTER :: f_pk(:) 
     13    TYPE(t_field), POINTER :: f_ps(:)  ! IN 
     14    TYPE(t_field), POINTER :: f_p(:)   ! IN 
     15    TYPE(t_field), POINTER :: f_pks(:) ! OUT 
     16    TYPE(t_field), POINTER :: f_pk(:)  ! OUT 
    1717   
    1818    REAL(rstd), POINTER :: ps(:) 
    … …  
    4949    INTEGER,INTENT(IN) :: offset 
    5050    INTEGER :: i,j,ij,l 
    51     REAL(rstd) :: alpha(iim*jjm,llm),beta(iim*jjm,llm) 
    52     REAL(rstd) :: delta  
    5351     
    54     IF(caldyn_exner == lmdz) THEN          ! LMD-Z style calculation of Exner pressure 
    55        !! Compute Alpha and Beta 
     52    ! surface : pks 
     53    IF (is_omp_level_master) THEN 
    5654        
    57        IF (is_omp_level_master) THEN 
    58        ! for llm layer 
    59          DO j=jj_begin-offset,jj_end+offset 
    60             DO i=ii_begin-offset,ii_end+offset 
    61                ij=(j-1)*iim+i 
    62                alpha(ij,llm) = 0. 
    63                beta (ij,llm) = 1./ (1+ 2*kappa) 
    64             ENDDO 
    65          ENDDO 
     55       DO j=jj_begin-offset,jj_end+offset 
     56          DO i=ii_begin-offset,ii_end+offset 
     57             ij=(j-1)*iim+i 
     58             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 
     59          ENDDO 
     60       ENDDO 
    6661        
    67        ! for other layer    
    68          DO l = llm-1 , 2 , -1 
    69             DO j=jj_begin-offset,jj_end+offset 
    70                DO i=ii_begin-offset,ii_end+offset 
    71                   ij=(j-1)*iim+i 
    72                   delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) ) 
    73                   alpha(ij,l)  = - p(ij,l+1) / delta * alpha(ij,l+1) 
    74                   beta (ij,l)  =   p(ij,l  ) / delta    
    75                ENDDO 
    76             ENDDO 
    77          ENDDO 
    78         
    79          !! Compute pk 
    80         
    81          ! for first layer 
    82          DO j=jj_begin-offset,jj_end+offset 
    83             DO i=ii_begin-offset,ii_end+offset 
    84                ij=(j-1)*iim+i 
    85                pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 
    86                pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )     /    & 
    87                     (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-(1.+2*kappa) )* p(ij,2) ) 
    88             ENDDO 
    89          ENDDO 
    90         
    91        ! for other layers 
    92         
    93          DO l = 2, llm 
    94             DO j=jj_begin-offset,jj_end+offset 
    95                DO i=ii_begin-offset,ii_end+offset 
    96                   ij=(j-1)*iim+i 
    97                   pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1) 
    98                ENDDO 
    99             ENDDO 
    100          ENDDO 
    101  
    102        ENDIF 
    103  
    104     ELSE ! Simple calculation of Exner pressure based on centered average 
    105        ! surface : pks 
    106        IF (is_omp_level_master) THEN 
    107  
    108          DO j=jj_begin-offset,jj_end+offset 
    109             DO i=ii_begin-offset,ii_end+offset 
    110                ij=(j-1)*iim+i 
    111                pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 
    112             ENDDO 
    113          ENDDO 
    114  
    115        ENDIF 
    116  
    117          ! 3D : pk 
    118        DO l = 1, llm 
    119           DO j=jj_begin-offset,jj_end+offset 
    120              DO i=ii_begin-offset,ii_end+offset 
    121                ij=(j-1)*iim+i 
    122                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa 
     62    ENDIF 
     63     
     64    ! 3D : pk 
     65    DO l = 1, llm 
     66       DO j=jj_begin-offset,jj_end+offset 
     67          DO i=ii_begin-offset,ii_end+offset 
     68             ij=(j-1)*iim+i 
     69             pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa 
    12370             ENDDO 
    12471          ENDDO 
    12572       ENDDO 
    126  
    127     END IF 
    128      
     73        
    12974  END SUBROUTINE compute_exner 
    13075   

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

    r354 r428  
    22  IMPLICIT NONE 
    33  PRIVATE 
    4  
     4   
    55  PUBLIC :: compute_geopotential 
    66CONTAINS 
     7   
     8  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_p,f_theta,f_phi) ! FIXME : never called, dry only 
     9    USE icosa 
     10    USE omp_para 
     11    USE pression_mod 
     12    USE theta2theta_rhodz_mod 
     13    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN                                        
     14         f_p(:), f_theta(:), f_phi(:)               ! OUT                                                          
     15    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:,:), theta_rhodz(:,:,:), & 
     16         phi(:,:), phis(:), ps(:) 
     17    INTEGER :: ind 
     18     
     19    DO ind=1,ndomain 
     20       IF (.NOT. assigned_domain(ind)) CYCLE 
     21       CALL swap_dimensions(ind) 
     22       CALL swap_geometry(ind) 
     23       ps = f_ps(ind) 
     24       p  = f_p(ind) 
     25!$OMP BARRIER                                                                                                      
     26       CALL compute_pression(ps,p,0) 
     27!$OMP BARRIER                                                                                                      
     28       theta_rhodz = f_theta_rhodz(ind) 
     29       theta = f_theta(ind) 
     30       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0) 
     31       phis = f_phis(ind) 
     32       phi = f_phi(ind) 
     33       CALL compute_geopotential(phis,p,theta,phi,0) 
     34    END DO 
    735 
    8   SUBROUTINE geopotential(f_phis,f_pks,f_pk,f_theta,f_phi) ! ORPHAN 
    9   USE icosa 
    10     TYPE(t_field), POINTER :: f_phis(:) 
    11     TYPE(t_field), POINTER :: f_pks(:) 
    12     TYPE(t_field), POINTER :: f_pk(:) 
    13     TYPE(t_field), POINTER :: f_theta(:) 
    14     TYPE(t_field), POINTER :: f_phi(:) 
    15    
    16     REAL(rstd), POINTER :: phis(:) 
    17     REAL(rstd), POINTER :: pks(:) 
    18     REAL(rstd), POINTER :: pk(:,:) 
    19     REAL(rstd), POINTER :: theta(:,:) 
    20     REAL(rstd), POINTER :: phi(:,:) 
    21     INTEGER :: ind 
     36  END SUBROUTINE thetarhodz2geopot 
    2237 
    23     DO ind=1,ndomain 
    24       IF (.NOT. assigned_domain(ind)) CYCLE 
    25       CALL swap_dimensions(ind) 
    26       CALL swap_geometry(ind) 
    27       phis=f_phis(ind) 
    28       pks=f_pks(ind) 
    29       pk=f_pk(ind) 
    30       theta=f_theta(ind) 
    31       phi=f_phi(ind) 
    32       CALL compute_geopotential(phis,pks,pk,theta,phi,0) 
    33     ENDDO 
    34    
    35   END SUBROUTINE geopotential 
    36    
    37   SUBROUTINE compute_geopotential(phis,pks,pk,theta,phi,offset) 
    38   USE icosa 
    39   USE omp_para 
    40   IMPLICIT NONE 
     38 
     39  SUBROUTINE compute_geopotential(phis,p,theta,phi,offset) 
     40    USE icosa 
     41    USE omp_para 
     42    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1) 
    4143    REAL(rstd),INTENT(IN) :: phis(iim*jjm) 
    42     REAL(rstd),INTENT(IN) :: pks(iim*jjm) 
    43     REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm) 
    44     REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm) 
     44    REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm,nqdyn) 
    4545    REAL(rstd),INTENT(OUT) :: phi(iim*jjm,llm) 
    4646    INTEGER,INTENT(IN) :: offset 
    4747    INTEGER :: i,j,ij,l 
     48    REAL(rstd) :: mg_ij, p_ij, exner_ij 
    4849 
    4950!!! Compute geopotential 
    … …  
    5758! flush pks, pk thetha, phis 
    5859!$OMP BARRIER 
    59    IF(is_omp_level_master) THEN 
    60      DO j=jj_begin-offset,jj_end+offset 
    61        DO i=ii_begin-offset,ii_end+offset 
    62          ij=(j-1)*iim+i 
    63          phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) ) 
    64        ENDDO 
    65      ENDDO 
     60    IF(is_omp_level_master) THEN 
     61       DO j=jj_begin-offset,jj_end+offset 
     62          DO i=ii_begin-offset,ii_end+offset 
     63             ij=(j-1)*iim+i 
     64             phi( ij,1 ) = phis( ij ) 
     65          END DO 
     66       END DO 
    6667           
    6768    ! for other layers 
    68      DO l = 2, llm 
     69     DO l = 1, llm-1 
    6970       DO j=jj_begin-offset,jj_end+offset 
    7071         DO i=ii_begin-offset,ii_end+offset 
    7172           ij=(j-1)*iim+i 
    72            phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  &  
    73                                          * (  pk(ij,l-1) -  pk(ij,l)    ) 
     73           mg_ij = p(ij,l+1)-p(ij,l) 
     74           ! v=RT/p=(kappa.cpp).Theta.(p/preff)**kappa /p 
     75           p_ij=.5*(p(ij,l+1)+p(ij,l)) 
     76           exner_ij = (p_ij/preff)**kappa ! NB : no cpp factor 
     77           phi(ij,l+1) = phi(ij,l) + (kappa*cpp)*mg_ij*theta(ij,l,1)*exner_ij/p_ij 
    7478         ENDDO 
    7579       ENDDO 

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

    r418 r428  
    158158  END SUBROUTINE write_output_fields_basic 
    159159   
    160   SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, & 
    161        f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p) 
    162     USE vorticity_mod 
    163     USE theta2theta_rhodz_mod 
    164     USE pression_mod 
    165     USE omega_mod 
    166     USE write_field_mod 
    167     USE vertical_interp_mod 
    168     USE wind_mod 
    169     TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), & 
    170          f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:) 
    171      
    172     REAL(rstd) :: out_pression_level 
    173     CHARACTER(LEN=255) :: str_pression 
    174     CHARACTER(LEN=255) :: physics_type 
    175      
    176     out_pression_level=0. 
    177     CALL getin("out_pression_level",out_pression_level)  
    178     WRITE(str_pression,*) INT(out_pression_level/100) 
    179     str_pression=ADJUSTL(str_pression) 
    180      
    181     CALL writefield("ps",f_ps) 
    182     CALL writefield("dps",f_dps) 
    183     CALL writefield("phis",f_phis) 
    184     CALL vorticity(f_u,f_buf_v) 
    185     CALL writefield("vort",f_buf_v) 
    186  
    187     CALL w_omega(f_ps, f_u, f_buf_i) 
    188     CALL writefield('omega', f_buf_i) 
    189     IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN 
    190       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level) 
    191       CALL writefield("omega"//TRIM(str_pression),f_buf_s) 
    192     ENDIF 
    193      
    194     ! Temperature 
    195 !    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME 
    196       
    197     CALL getin('physics',physics_type) 
    198     IF (TRIM(physics_type)=='dcmip') THEN 
    199       CALL Tv2T(f_buf_i,f_q,f_buf1_i)  
    200       CALL writefield("T",f_buf1_i) 
    201       IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN 
    202         CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level) 
    203         CALL writefield("T"//TRIM(str_pression),f_buf_s) 
    204       ENDIF 
    205     ELSE 
    206       CALL writefield("T",f_buf_i) 
    207       IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN 
    208         CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level) 
    209         CALL writefield("T"//TRIM(str_pression),f_buf_s) 
    210       ENDIF 
    211     ENDIF 
    212     
    213     ! velocity components 
    214     CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i) 
    215     CALL writefield("ulon",f_buf1_i) 
    216     CALL writefield("ulat",f_buf2_i) 
    217  
    218     IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN 
    219       CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level) 
    220       CALL writefield("ulon"//TRIM(str_pression),f_buf_s) 
    221       CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level) 
    222       CALL writefield("ulat"//TRIM(str_pression),f_buf_s) 
    223     ENDIF 
    224      
    225     ! geopotential ! FIXME 
    226     CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i) 
    227     CALL writefield("p",f_buf_p) 
    228 !    CALL writefield("phi",f_geopot)   ! geopotential 
    229     CALL writefield("theta",f_buf1_i) ! potential temperature 
    230     CALL writefield("pk",f_buf2_i)    ! Exner pressure 
    231    
    232   END SUBROUTINE write_output_fields 
    233  
    234 !------------------- Conversion from prognostic to observable variables ------------------ 
     160 !------------------- Conversion from prognostic to observable variables ------------------ 
    235161 
    236162  SUBROUTINE progonostic_vel_to_horiz(f_geopot, f_ps, f_rhodz, f_u, f_W, f_uh, f_uz) 
    … …  
    308234  END SUBROUTINE compute_prognostic_vel_to_horiz 
    309235 
    310   SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi)  
    311     USE field_mod 
    312     USE pression_mod 
    313     USE exner_mod 
    314     USE geopotential_mod 
    315     USE theta2theta_rhodz_mod 
    316     TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN 
    317          f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT 
    318     REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:,:), & 
    319          phi(:,:), phis(:), ps(:), pks(:) 
    320     INTEGER :: ind 
    321      
    322     DO ind=1,ndomain 
    323        IF (.NOT. assigned_domain(ind)) CYCLE 
    324        CALL swap_dimensions(ind) 
    325        CALL swap_geometry(ind) 
    326        ps = f_ps(ind) 
    327        p  = f_p(ind) 
    328 !$OMP BARRIER 
    329        CALL compute_pression(ps,p,0) 
    330        pk = f_pk(ind) 
    331        pks = f_pks(ind) 
    332 !$OMP BARRIER 
    333        CALL compute_exner(ps,p,pks,pk,0) 
    334 !$OMP BARRIER 
    335        theta_rhodz = f_theta_rhodz(ind) 
    336        theta = f_theta(ind) 
    337        CALL compute_theta_rhodz2theta(ps, theta_rhodz(:,:,1),theta,0) 
    338        phis = f_phis(ind) 
    339        phi = f_phi(ind) 
    340        CALL compute_geopotential(phis,pks,pk,theta,phi,0) 
    341     END DO 
    342  
    343   END SUBROUTINE thetarhodz2geopot 
    344    
    345236  SUBROUTINE Tv2T(f_Tv, f_q, f_T) 
    346237    TYPE(t_field), POINTER :: f_TV(:) 

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

    r387 r428  
    11MODULE theta2theta_rhodz_mod 
    22  USE field_mod 
    3    
     3  PRIVATE 
    44  TYPE(t_field), POINTER, SAVE  :: f_p(:) 
    5   TYPE(t_field), POINTER, SAVE  :: f_pks(:) 
    6   TYPE(t_field), POINTER, SAVE  :: f_pk(:) 
    7  
    8   PRIVATE :: f_p,f_pk,f_pks  
     5 
     6  PUBLIC :: init_theta2theta_rhodz, theta_rhodz2theta, & 
     7       theta_rhodz2temperature, temperature2theta_rhodz, & 
     8       theta2theta_rhodz, & 
     9       compute_theta2theta_rhodz, compute_theta_rhodz2theta 
    910 
    1011CONTAINS 
    … …  
    1415  USE field_mod 
    1516  IMPLICIT NONE 
    16     CALL allocate_field(f_p,field_t,type_real,llm+1,name='p (theta2theta_rhodz_mod)') 
    17     CALL allocate_field(f_pk,field_t,type_real,llm,name='pk (theta2theta_rhodz_mod)') 
    18     CALL allocate_field(f_pks,field_t,type_real,name='pks (theta2theta_rhodz_mod)') 
    19      
     17    CALL allocate_field(f_p,field_t,type_real,llm+1,name='p (theta2theta_rhodz_mod)')     
    2018  END SUBROUTINE init_theta2theta_rhodz 
    21  
    2219 
    2320 
    … …  
    6158    REAL(rstd), POINTER :: temp(:,:) 
    6259    REAL(rstd), POINTER :: p(:,:) 
    63     REAL(rstd), POINTER :: pk(:,:) 
    64     REAL(rstd), POINTER :: pks(:) 
    6560    INTEGER :: ind 
    6661 
    … …  
    7166      ps=f_ps(ind) 
    7267      p=f_p(ind) 
    73       pks=f_pks(ind) 
    74       pk=f_pk(ind) 
    7568      theta_rhodz=f_theta_rhodz(ind) 
    7669      temp=f_temp(ind) 
    … …  
    7972      CALL compute_pression(ps,p,0) 
    8073!$OMP BARRIER 
    81       CALL compute_exner(ps,p,pks,pk,0) 
    82 !$OMP BARRIER 
    83       CALL compute_theta_rhodz2temperature(p, pk, theta_rhodz(:,:,1),temp,0) 
     74      CALL compute_theta_rhodz2temperature(p, theta_rhodz(:,:,1),temp,0) 
    8475    ENDDO 
    8576!$OMP BARRIER 
    … …  
    10091    REAL(rstd), POINTER :: temp(:,:) 
    10192    REAL(rstd), POINTER :: p(:,:) 
    102     REAL(rstd), POINTER :: pk(:,:) 
    103     REAL(rstd), POINTER :: pks(:) 
    10493    INTEGER :: ind 
    10594 
    … …  
    11099      ps=f_ps(ind) 
    111100      p=f_p(ind) 
    112       pks=f_pks(ind) 
    113       pk=f_pk(ind) 
    114101      theta_rhodz=f_theta_rhodz(ind) 
    115102      temp=f_temp(ind) 
    … …  
    118105      CALL compute_pression(ps,p,0) 
    119106!$OMP BARRIER 
    120       CALL compute_exner(ps,p,pks,pk,0) 
    121 !$OMP BARRIER 
    122       CALL compute_temperature2theta_rhodz(p, pk, temp, theta_rhodz(:,:,1), 0) 
     107      CALL compute_temperature2theta_rhodz(p, temp, theta_rhodz(:,:,1), 0) 
    123108    ENDDO 
    124109!$OMP BARRIER 
    … …  
    213198 
    214199 
    215   SUBROUTINE compute_theta_rhodz2temperature(p,pk,theta_rhodz,temp,offset) 
     200  SUBROUTINE compute_theta_rhodz2temperature(p,theta_rhodz,temp,offset) 
    216201  USE icosa 
    217202  USE pression_mod 
    … …  
    220205  IMPLICIT NONE 
    221206    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1) 
    222     REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm) 
    223207    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm) 
    224208    REAL(rstd),INTENT(OUT) :: temp(iim*jjm,llm) 
    225209    INTEGER,INTENT(IN) :: offset 
     210    REAL(rstd) :: pk_ij 
    226211    INTEGER :: i,j,ij,l 
    227212         
    … …  
    232217        DO i=ii_begin-offset,ii_end+offset 
    233218          ij=(j-1)*iim+i 
    234           temp(ij,l) = ( theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g) ) * pk(ij,l) / cpp 
     219          pk_ij=((.5/preff)*(p(ij,l)+p(ij,l+1)))**kappa 
     220          temp(ij,l) = ( theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g) ) * pk_ij 
    235221        ENDDO 
    236222      ENDDO 
    … …  
    241227  END SUBROUTINE compute_theta_rhodz2temperature 
    242228 
    243   SUBROUTINE compute_temperature2theta_rhodz(p,pk,temp,theta_rhodz,offset) 
     229  SUBROUTINE compute_temperature2theta_rhodz(p,temp,theta_rhodz,offset) 
    244230  USE icosa 
    245231  USE pression_mod 
    … …  
    248234  IMPLICIT NONE 
    249235    REAL(rstd),INTENT(IN)  :: p(iim*jjm,llm+1) 
    250     REAL(rstd),INTENT(IN)  :: pk(iim*jjm,llm) 
    251236    REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 
    252237    REAL(rstd),INTENT(IN)  :: temp(iim*jjm,llm) 
    253238    INTEGER,INTENT(IN) :: offset 
     239    REAL(rstd) :: pk_ij 
    254240    INTEGER :: i,j,ij,l 
    255241 
    … …  
    262248        DO i=ii_begin-offset,ii_end+offset 
    263249          ij=(j-1)*iim+i 
    264           theta_rhodz(ij,l) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp ) 
     250          pk_ij=((.5/preff)*(p(ij,l)+p(ij,l+1)))**kappa 
     251          theta_rhodz(ij,l) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / pk_ij 
    265252        ENDDO 
    266253      ENDDO 
Note: See TracChangeset for help on using the changeset viewer.