Changeset 405


Ignore:
Timestamp:
06/08/16 22:24:19 (8 years ago)
Author:
dubos
Message:

Moist dynamics (to be tested)

Location:
codes/icosagcm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • codes/icosagcm/trunk/arch/arch-X64_YELLOWSTONE.fcm

    r399 r405  
    77%FPP_DEF             KEY_NONE CPP_NETCDF4 CPP_NF90_DEF_VAR_CHUNKING 
    88%BASE_FFLAGS         -i4 -r8 -auto -align all -I${MKLROOT}/include 
    9 %PROD_FFLAGS         -g -traceback -O3 -vec-report2 
     9%PROD_FFLAGS         -g -traceback -O3 
    1010%DEV_FFLAGS          -g -O1 -traceback 
    1111%DEBUG_FFLAGS        -g -traceback -check bounds -fp-model strict 
  • codes/icosagcm/trunk/src/caldyn_kernels_base.f90

    r404 r405  
    3737 
    3838    INTEGER :: i,j,ij,l 
    39     REAL(rstd) :: Rd, p_ik, exner_ik, temp_ik 
     39    REAL(rstd) :: Rd, p_ik, exner_ik, temp_ik, qv, chi, Rmix 
    4040    INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext 
    4141 
     
    7979       ! uppermost layer 
    8080        
    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         
     81       SELECT CASE(caldyn_thermo) 
     82          CASE(thermo_theta, thermo_entropy) 
     83             !DIR$ SIMD 
     84             DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     85                pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 
     86             END DO 
     87             ! other layers 
     88             DO l = llm-1, 1, -1 
     89                !DIR$ SIMD 
     90                DO ij=ij_omp_begin_ext,ij_omp_end_ext          
     91                   pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1)) 
     92                END DO 
     93             END DO 
     94             ! surface pressure (for diagnostics) 
     95             IF(caldyn_eta==eta_lag) THEN 
     96                DO ij=ij_omp_begin_ext,ij_omp_end_ext          
     97                   ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 
     98                END DO 
     99             END IF 
     100          CASE(thermo_moist) ! theta(ij,l,2) = qv = mv/md 
     101             !DIR$ SIMD 
     102             DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     103                pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)*(1.+theta(ij,l,2)) 
     104             END DO 
     105             ! other layers 
     106             DO l = llm-1, 1, -1 
     107                !DIR$ SIMD 
     108                DO ij=ij_omp_begin_ext,ij_omp_end_ext          
     109                   pk(ij,l) = pk(ij,l+1) + (.5*g)*(          & 
     110                        rhodz(ij,l)  *(1.+theta(ij,l,2)) +   & 
     111                        rhodz(ij,l+1)*(1.+theta(ij,l+1,2)) ) 
     112                END DO 
     113             END DO 
     114             ! surface pressure (for diagnostics) 
     115             IF(caldyn_eta==eta_lag) THEN 
     116                DO ij=ij_omp_begin_ext,ij_omp_end_ext          
     117                   ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)*(1.+theta(ij,l,2)) 
     118                END DO 
     119             END IF 
     120          END SELECT 
     121 
    99122       DO l = 1,llm 
    100123          SELECT CASE(caldyn_thermo) 
     
    117140                geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik 
    118141             ENDDO 
     142          CASE(thermo_moist) ! theta is moist pseudo-entropy per dry air mass 
     143             DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     144                p_ik = pk(ij,l) 
     145                qv = theta(ij,l,2) ! water vaper mixing ratio = mv/md 
     146                Rmix = Rd+qv*Rv 
     147                chi = ( theta(ij,l,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff) 
     148                temp_ik = Treff*exp(chi) 
     149                pk(ij,l) = temp_ik 
     150                ! specific volume v = R*T/p = dphi/g/rhodz 
     151                ! R = (Rd + qv.Rv)/(1+qv) 
     152                geopot(ij,l+1) = geopot(ij,l) + g*Rmix*rhodz(ij,l)*temp_ik/(p_ik*(1+qv)) 
     153             ENDDO 
    119154          CASE DEFAULT 
    120155             STOP 
     
    132167  SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du) 
    133168    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
    134     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm) 
     169    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn) 
    135170    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    136171    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence 
     
    139174    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1) 
    140175    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm) 
    141     REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm) 
     176    REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn) 
    142177    REAL(rstd),INTENT(OUT) :: dps(iim*jjm) 
    143178 
    144179    ! temporary variable     
    145     INTEGER :: i,j,ij,l 
     180    INTEGER :: i,j,ij,l,iq 
    146181    REAL(rstd) :: p_ik, exner_ik 
    147182    INTEGER    :: ij_omp_begin, ij_omp_end 
     
    204239    !$OMP BARRIER 
    205240 
    206     DO l=ll_begin,ll_endm1 
    207        !DIR$ SIMD 
    208        DO ij=ij_begin,ij_end          
    209           dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1))) 
    210        ENDDO 
    211     ENDDO 
    212  
    213     DO l=ll_beginp1,ll_end 
    214        !DIR$ SIMD 
    215        DO ij=ij_begin,ij_end          
    216           dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) ) 
    217        ENDDO 
    218     ENDDO 
    219  
     241    DO iq=1,nqdyn 
     242       DO l=ll_begin,ll_endm1 
     243       !DIR$ SIMD 
     244          DO ij=ij_begin,ij_end          
     245             dtheta_rhodz(ij, l, iq) = dtheta_rhodz(ij, l, iq)  -   0.5 * & 
     246                  ( wflux(ij,l+1) * (theta(ij,l,iq) + theta(ij,l+1,iq))) 
     247          END DO 
     248       END DO 
     249       DO l=ll_beginp1,ll_end 
     250          !DIR$ SIMD 
     251          DO ij=ij_begin,ij_end          
     252             dtheta_rhodz(ij, l, iq) = dtheta_rhodz(ij, l, iq)  +   0.5 * & 
     253                  ( wflux(ij,l) * (theta(ij,l-1,iq) + theta(ij,l,iq) ) ) 
     254          END DO 
     255       END DO 
     256    END DO 
    220257 
    221258    ! Compute vertical transport 
  • codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90

    r404 r405  
    382382    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! OUT if tau>0 
    383383    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm) 
    384     REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm) 
     384    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) 
    385385    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) 
    386386    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) 
    387387    REAL(rstd),INTENT(INOUT)   :: du(iim*3*jjm,llm) 
    388388    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
     389    REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function 
    389390 
    390391    INTEGER :: i,j,ij,l 
    391     REAL(rstd) :: due_right, due_lup, due_ldown 
     392    REAL(rstd) :: Rd, qv, temp, chi, nu, due_right, due_lup, due_ldown 
    392393 
    393394    CALL trace_start("compute_caldyn_fast") 
     395 
     396    Rd=cpp*kappa 
    394397 
    395398    ! Compute Bernoulli term 
    396399    IF(boussinesq) THEN 
    397  
    398400       DO l=ll_begin,ll_end 
    399401          !DIR$ SIMD 
     
    417419             DO ij=ij_begin,ij_end          
    418420                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
    419                      + pk(ij,l)*(cpp-theta(ij,l)) ! pk=temperature, theta=entropy 
     421                     + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy 
    420422             END DO 
     423          CASE(thermo_moist)  
     424             !DIR$ SIMD 
     425             DO ij=ij_begin,ij_end          
     426                ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T) 
     427                ! Bd = Phi + gibbs_d 
     428                ! Bv = Phi + gibbs_v 
     429                ! pk=temperature, theta=entropy 
     430                qv = theta(ij,l,2) 
     431                temp = pk(ij,l) 
     432                chi = log(temp/Treff) 
     433                nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff) 
     434                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
     435                     + temp*(cpp*(1.-chi)+Rd*nu) 
     436                berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
     437                     + temp*(cppv*(1.-chi)+Rv*nu) 
     438            END DO 
    421439          END SELECT 
    422440       END DO 
     
    426444!!! u:=u+tau*du, du = -grad(B)-theta.grad(pi) 
    427445    DO l=ll_begin,ll_end 
    428        !DIR$ SIMD 
    429        DO ij=ij_begin,ij_end          
    430           due_right =  0.5*(theta(ij,l)+theta(ij+t_right,l))  & 
    431                           *(pk(ij+t_right,l)-pk(ij,l))        & 
    432                          +  berni(ij+t_right,l)-berni(ij,l) 
    433           due_lup = 0.5*(theta(ij,l)+theta(ij+t_lup,l))    & 
    434                        *(pk(ij+t_lup,l)-pk(ij,l))          & 
    435                          +  berni(ij+t_lup,l)-berni(ij,l) 
    436           due_ldown = 0.5*(theta(ij,l)+theta(ij+t_ldown,l)) & 
    437                          *(pk(ij+t_ldown,l)-pk(ij,l))      & 
    438                         +   berni(ij+t_ldown,l)-berni(ij,l) 
    439           du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right 
    440           du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup 
    441           du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown 
    442           u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l) 
    443           u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l) 
    444           u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l) 
    445        ENDDO 
    446     ENDDO 
    447  
     446       IF(caldyn_thermo == thermo_moist) THEN 
     447          !DIR$ SIMD 
     448          DO ij=ij_begin,ij_end          
     449             due_right =  berni(ij+t_right,l)-berni(ij,l)       & 
     450                  + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   & 
     451                       *(pk(ij+t_right,l)-pk(ij,l))             & 
     452                  + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   & 
     453                       *(berniv(ij+t_right,l)-berniv(ij,l)) 
     454              
     455             due_lup = berni(ij+t_lup,l)-berni(ij,l)            & 
     456                  + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     & 
     457                       *(pk(ij+t_lup,l)-pk(ij,l))               & 
     458                  + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     & 
     459                       *(berniv(ij+t_lup,l)-berniv(ij,l)) 
     460              
     461             due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        & 
     462                  + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   & 
     463                       *(pk(ij+t_ldown,l)-pk(ij,l))             & 
     464                  + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   & 
     465                       *(berniv(ij+t_ldown,l)-berniv(ij,l)) 
     466              
     467             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right 
     468             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup 
     469             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown 
     470             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l) 
     471             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l) 
     472             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l) 
     473          END DO 
     474       ELSE 
     475          !DIR$ SIMD 
     476          DO ij=ij_begin,ij_end          
     477             due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  & 
     478                  *(pk(ij+t_right,l)-pk(ij,l))        & 
     479                  +  berni(ij+t_right,l)-berni(ij,l) 
     480             due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    & 
     481                  *(pk(ij+t_lup,l)-pk(ij,l))          & 
     482                  +  berni(ij+t_lup,l)-berni(ij,l) 
     483             due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) & 
     484                  *(pk(ij+t_ldown,l)-pk(ij,l))      & 
     485                  +   berni(ij+t_ldown,l)-berni(ij,l) 
     486             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right 
     487             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup 
     488             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown 
     489             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l) 
     490             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l) 
     491             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l) 
     492          END DO 
     493       END IF 
     494    END DO 
     495        
    448496    CALL trace_end("compute_caldyn_fast") 
    449497 
Note: See TracChangeset for help on using the changeset viewer.