Changeset 404


Ignore:
Timestamp:
06/08/16 18:18:08 (8 years ago)
Author:
dubos
Message:

Towards moist dynamics - tested with DCMIP41

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

Legend:

Unmodified
Added
Removed

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

    r387 r404  
    5454 
    5555! temporary shared variable 
    56     REAL(rstd),POINTER  :: theta(:,:)   
     56    REAL(rstd),POINTER  :: theta(:,:,:)   
    5757    REAL(rstd),POINTER  :: pk(:,:) 
    5858    REAL(rstd),POINTER  :: geopot(:,:) 
    … …  
    107107       mass=f_mass(ind) 
    108108       theta = f_theta(ind) 
    109        CALL compute_theta(ps,theta_rhodz(:,:,1), mass,theta) 
     109       CALL compute_theta(ps,theta_rhodz, mass,theta) 
    110110       pk = f_pk(ind) 
    111111       geopot = f_geopot(ind) 
    … …  
    113113       IF(hydrostatic) THEN 
    114114          du(:,:)=0. 
    115           CALL compute_geopot(ps,mass,theta, pk,geopot) 
     115          CALL compute_geopot(mass,theta, ps,pk,geopot) 
    116116       ELSE 
    117117          W = f_W(ind) 
    … …  
    163163          CALL compute_caldyn_slow_NH(u,mass,geopot,W, hflux,du,dPhi,dW) 
    164164       END IF 
    165        CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz(:,:,1),du) 
     165       CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du) 
    166166       IF(caldyn_eta==eta_mass) THEN 
    167167          wflux=f_wflux(ind) 

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

    r401 r404  
    2929  !**************************** Geopotential ***************************** 
    3030 
    31   SUBROUTINE compute_geopot(ps,rhodz,theta, pk,geopot)  
     31  SUBROUTINE compute_geopot(rhodz,theta, ps,pk,geopot)  
     32    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm) 
     33    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) ! active scalars : theta/entropy, moisture, ... 
    3234    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm) 
    33     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm) 
    34     REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)    ! potential temperature 
    3535    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)       ! Exner function (compressible) /Lagrange multiplier (Boussinesq) 
    3636    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential 
    … …  
    4848    Rd = kappa*cpp 
    4949 
    50     IF(caldyn_eta==eta_mass .AND. .NOT. DEC) THEN 
    51 !!! Compute exner function and geopotential 
    52        STOP 
     50    ! Pressure is computed first top-down (temporarily stored in pk) 
     51    ! Then Exner pressure and geopotential are computed bottom-up 
     52    ! Works also when caldyn_eta=eta_mass           
     53 
     54    IF(boussinesq) THEN ! compute geopotential and pk=Lagrange multiplier 
     55       ! specific volume 1 = dphi/g/rhodz 
     56       !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency 
    5357       DO l = 1,llm 
    5458          !DIR$ SIMD 
    55           DO ij=ij_omp_begin_ext,ij_omp_end_ext 
    56              p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 
    57              exner_ik = cpp * (p_ik/preff) ** kappa 
    58              pk(ij,l) = exner_ik 
    59              ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
    60              geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
     59          DO ij=ij_omp_begin_ext,ij_omp_end_ext          
     60             geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l) 
    6161          ENDDO 
    6262       ENDDO 
    63     ELSE  
    64        ! We are using DEC or a Lagrangian vertical coordinate 
    65        ! Pressure is computed first top-down (temporarily stored in pk) 
    66        ! Then Exner pressure and geopotential are computed bottom-up 
    67        ! Works also when caldyn_eta=eta_mass           
    68  
    69        IF(boussinesq) THEN ! compute geopotential and pk=Lagrange multiplier 
    70           ! specific volume 1 = dphi/g/rhodz 
    71           !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency 
    72           DO l = 1,llm 
    73              !DIR$ SIMD 
    74              DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    75                 geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l) 
    76              ENDDO 
    77           ENDDO 
    78           ! use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure)  
    79           ! uppermost layer 
     63       ! use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure)  
     64       ! uppermost layer 
     65       !DIR$ SIMD 
     66       DO ij=ij_begin_ext,ij_end_ext          
     67          pk(ij,llm) = ptop + (.5*g)*theta(ij,llm,1)*rhodz(ij,llm) 
     68       END DO 
     69       ! other layers 
     70       DO l = llm-1, 1, -1 
     71          !          !$OMP DO SCHEDULE(STATIC)  
    8072          !DIR$ SIMD 
    8173          DO ij=ij_begin_ext,ij_end_ext          
    82              pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm) 
     74             pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l,1)*rhodz(ij,l)+theta(ij,l+1,1)*rhodz(ij,l+1)) 
    8375          END DO 
    84           ! other layers 
    85           DO l = llm-1, 1, -1 
    86              !          !$OMP DO SCHEDULE(STATIC)  
     76       END DO 
     77       ! now pk contains the Lagrange multiplier (pressure) 
     78    ELSE ! non-Boussinesq, compute pressure, Exner pressure or temperature, then geopotential 
     79       ! uppermost layer 
     80        
     81       !DIR$ SIMD 
     82       DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     83          pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 
     84       END DO 
     85       ! other layers 
     86       DO l = llm-1, 1, -1 
     87          !DIR$ SIMD 
     88          DO ij=ij_omp_begin_ext,ij_omp_end_ext          
     89             pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1)) 
     90          END DO 
     91       END DO 
     92       ! surface pressure (for diagnostics) 
     93       IF(caldyn_eta==eta_lag) THEN 
     94          DO ij=ij_omp_begin_ext,ij_omp_end_ext          
     95             ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 
     96          END DO 
     97       END IF 
     98        
     99       DO l = 1,llm 
     100          SELECT CASE(caldyn_thermo) 
     101          CASE(thermo_theta) 
    87102             !DIR$ SIMD 
    88              DO ij=ij_begin_ext,ij_end_ext          
    89                 pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1)) 
    90              END DO 
    91           END DO 
    92           ! now pk contains the Lagrange multiplier (pressure) 
    93        ELSE ! non-Boussinesq, compute pressure, Exner pressure or temperature, then geopotential 
    94           ! uppermost layer 
    95  
    96           !DIR$ SIMD 
    97           DO ij=ij_omp_begin_ext,ij_omp_end_ext 
    98              pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 
    99           END DO 
    100           ! other layers 
    101           DO l = llm-1, 1, -1 
     103             DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     104                p_ik = pk(ij,l) 
     105                exner_ik = cpp * (p_ik/preff) ** kappa 
     106                pk(ij,l) = exner_ik 
     107                ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
     108                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l,1)*exner_ik/p_ik 
     109             ENDDO 
     110          CASE(thermo_entropy) ! theta is in fact entropy = cpp*log(theta/Treff) = cpp*log(T/Treff) - Rd*log(p/preff) 
    102111             !DIR$ SIMD 
    103              DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    104                 pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1)) 
    105              END DO 
    106           END DO 
    107           ! surface pressure (for diagnostics) 
    108           IF(caldyn_eta==eta_lag) THEN 
    109              DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    110                 ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 
    111              END DO 
    112           END IF 
    113  
    114           DO l = 1,llm 
    115              SELECT CASE(caldyn_thermo) 
    116              CASE(thermo_theta) 
    117                 !DIR$ SIMD 
    118                 DO ij=ij_omp_begin_ext,ij_omp_end_ext 
    119                    p_ik = pk(ij,l) 
    120                    exner_ik = cpp * (p_ik/preff) ** kappa 
    121                    pk(ij,l) = exner_ik 
    122                    ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
    123                    geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
    124                 ENDDO 
    125              CASE(thermo_entropy) ! theta is in fact entropy = cpp*log(theta/Treff) = cpp*log(T/Treff) - Rd*log(p/preff) 
    126                 !DIR$ SIMD 
    127                 DO ij=ij_omp_begin_ext,ij_omp_end_ext 
    128                    p_ik = pk(ij,l) 
    129                    temp_ik = Treff*exp((theta(ij,l) + Rd*log(p_ik/preff))/cpp) 
    130                    pk(ij,l) = temp_ik 
    131                    ! specific volume v = Rd*T/p = dphi/g/rhodz 
    132                    geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik 
    133                 ENDDO 
    134              CASE DEFAULT 
    135                 STOP 
    136              END SELECT 
    137           ENDDO 
    138        END IF 
    139  
     112             DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     113                p_ik = pk(ij,l) 
     114                temp_ik = Treff*exp((theta(ij,l,1) + Rd*log(p_ik/preff))/cpp) 
     115                pk(ij,l) = temp_ik 
     116                ! specific volume v = Rd*T/p = dphi/g/rhodz 
     117                geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik 
     118             ENDDO 
     119          CASE DEFAULT 
     120             STOP 
     121          END SELECT 
     122       ENDDO 
    140123    END IF 
    141124 

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

    r401 r404  
    2020 
    2121  SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta) 
    22     REAL(rstd),INTENT(IN)  :: ps(iim*jjm) 
    23     REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm) 
     22    REAL(rstd),INTENT(IN)    :: ps(iim*jjm) 
     23    REAL(rstd),INTENT(IN)    :: theta_rhodz(iim*jjm,llm,nqdyn) 
    2424    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm) 
    25     REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm) 
    26     INTEGER :: ij,l 
     25    REAL(rstd),INTENT(OUT)   :: theta(iim*jjm,llm,nqdyn) 
     26    INTEGER :: ij,l,iq 
    2727    REAL(rstd) :: m 
    28     CALL trace_start("compute_theta")   
    29  
    30     IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta 
     28    CALL trace_start("compute_theta") 
     29 
     30    IF(caldyn_eta==eta_mass) THEN ! Compute mass 
    3131       DO l = ll_begin,ll_end 
    3232          !DIR$ SIMD 
    … …  
    3838             END IF 
    3939             rhodz(ij,l) = m/g 
    40              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    41           ENDDO 
    42        ENDDO 
    43     ELSE ! Compute only theta 
    44        DO l = ll_begin,ll_end 
     40          END DO 
     41       END DO 
     42    END IF 
     43 
     44    DO l = ll_begin,ll_end 
     45       DO iq=1,nqdyn 
    4546          !DIR$ SIMD 
    4647          DO ij=ij_begin_ext,ij_end_ext 
    47              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    48           ENDDO 
    49        ENDDO 
    50     END IF 
     48             theta(ij,l,iq) = theta_rhodz(ij,l,iq)/rhodz(ij,l) 
     49          END DO 
     50       END DO 
     51    END DO 
    5152 
    5253    CALL trace_end("compute_theta") 
    … …  
    451452  SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du) 
    452453    REAL(rstd),INTENT(IN)  :: hflux(3*iim*jjm,llm)  ! hflux in kg/s 
    453     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature 
     454    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn) ! active scalars 
    454455    REAL(rstd),INTENT(IN)  :: qu(3*iim*jjm,llm) 
    455456    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence 
    456     REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) 
     457    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn) 
    457458    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm) 
    458459     
    459460    REAL(rstd) :: Ftheta(3*iim*jjm)  ! potential temperature flux 
    460461    REAL(rstd) :: uu_right, uu_lup, uu_ldown 
    461     INTEGER :: ij,l,kdown 
     462    INTEGER :: ij,iq,l,kdown 
    462463 
    463464    CALL trace_start("compute_caldyn_Coriolis") 
    … …  
    466467       ! compute theta flux 
    467468       !DIR$ SIMD 
    468        DO ij=ij_begin_ext,ij_end_ext       
    469           Ftheta(ij+u_right) = 0.5*(theta(ij,l)+theta(ij+t_right,l)) & 
     469       DO iq=1,nqdyn 
     470          DO ij=ij_begin_ext,ij_end_ext       
     471             Ftheta(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) & 
    470472                                  * hflux(ij+u_right,l) 
    471           Ftheta(ij+u_lup)   = 0.5*(theta(ij,l)+theta(ij+t_lup,l)) & 
    472                                   * hflux(ij+u_lup,l) 
    473           Ftheta(ij+u_ldown) = 0.5*(theta(ij,l)+theta(ij+t_ldown,l)) & 
    474                                   * hflux(ij+u_ldown,l) 
    475        ENDDO 
    476        ! compute horizontal divergence of fluxes 
     473             Ftheta(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) & 
     474                  * hflux(ij+u_lup,l) 
     475             Ftheta(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) & 
     476                  * hflux(ij+u_ldown,l) 
     477          END DO 
     478          ! horizontal divergence of fluxes 
     479          DO ij=ij_begin,ij_end          
     480             ! dtheta_rhodz =  -div(flux.theta) 
     481             dtheta_rhodz(ij,l,iq)= & 
     482                  -1./Ai(ij)*(ne_right*Ftheta(ij+u_right)    +  & 
     483                  ne_rup*Ftheta(ij+u_rup)        +  &   
     484                  ne_lup*Ftheta(ij+u_lup)        +  &   
     485                  ne_left*Ftheta(ij+u_left)      +  &   
     486                  ne_ldown*Ftheta(ij+u_ldown)    +  & 
     487                  ne_rdown*Ftheta(ij+u_rdown) ) 
     488          END DO 
     489       END DO 
     490 
    477491       DO ij=ij_begin,ij_end          
    478492          ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 
    … …  
    482496               ne_left*hflux(ij+u_left,l)     +  &   
    483497               ne_ldown*hflux(ij+u_ldown,l)   +  &   
    484                ne_rdown*hflux(ij+u_rdown,l))     
    485           ! dtheta_rhodz =  -div(flux.theta) 
    486           dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right)    +  & 
    487                ne_rup*Ftheta(ij+u_rup)        +  &   
    488                ne_lup*Ftheta(ij+u_lup)        +  &   
    489                ne_left*Ftheta(ij+u_left)      +  &   
    490                ne_ldown*Ftheta(ij+u_ldown)    +  &   
    491                ne_rdown*Ftheta(ij+u_rdown)) 
    492        ENDDO 
    493     END DO 
     498               ne_rdown*hflux(ij+u_rdown,l)) 
     499       END DO ! ij 
     500    END DO ! llm 
    494501 
    495502!!! Compute potential vorticity (Coriolis) contribution to du 

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

    r401 r404  
    99  REAL(rstd),SAVE :: kappa=0.2857143 
    1010  REAL(rstd),SAVE :: cpp=1004.70885 
    11   REAL(rstd),SAVE :: cppv=1004.70885 ! FIXME 
     11  REAL(rstd),SAVE :: cppv=1860. 
    1212  REAL(rstd),SAVE :: Rv=461.5 
    1313  REAL(rstd),SAVE :: Treff=273. 
    … …  
    3939    CALL getin("kappa",kappa)   
    4040    CALL getin("cpp",cpp)   
     41    CALL getin("cppv",cppv) 
     42    CALL getin("Rv",Rv) 
    4143    CALL getin("preff",preff)   
    4244    CALL getin("Treff",Treff)   

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

    r403 r404  
    3333    CALL allocate_field(f_buf_s,   field_t,type_real, name="buf_s") 
    3434 
    35     CALL allocate_field(f_theta, field_t,type_real,llm,  name='theta')  ! potential temperature 
    36     CALL allocate_field(f_pmid,  field_t,type_real,llm,  name='pmid')   ! mid layer pressure 
     35    CALL allocate_field(f_theta, field_t,type_real,llm,nqdyn,  name='theta') ! potential temperature 
     36    CALL allocate_field(f_pmid,  field_t,type_real,llm,  name='pmid')       ! mid layer pressure 
    3737  END SUBROUTINE init_observable 
    3838   
Note: See TracChangeset for help on using the changeset viewer.