Changeset 854 for codes/icosagcm/devel


Ignore:
Timestamp:
05/05/19 21:39:45 (5 years ago)
Author:
jisesh
Message:

devel: separate module for compute_caldyn_fast and removed nu duplication in it

Location:
codes/icosagcm/devel/src/dynamics
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed
  • codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90

    r853 r854  
    1111  USE compute_caldyn_slow_NH_mod, ONLY : compute_caldyn_slow_NH 
    1212  USE compute_caldyn_solver_mod, ONLY : compute_caldyn_solver 
     13  USE compute_caldyn_fast_mod, ONLY : compute_caldyn_fast 
    1314  IMPLICIT NONE 
    1415  PRIVATE 
  • codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90

    r853 r854  
    1313  LOGICAL, SAVE :: debug_hevi_solver = .FALSE. 
    1414 
    15   PUBLIC :: compute_caldyn_fast,compute_NH_geopot 
     15  PUBLIC :: compute_NH_geopot 
    1616 
    1717CONTAINS 
     
    183183  END SUBROUTINE compute_NH_geopot 
    184184 
    185   SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du) 
    186     REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs 
    187     REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! OUT if tau>0 
    188     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm) 
    189     REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) 
    190     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) 
    191     REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) 
    192     REAL(rstd),INTENT(INOUT)   :: du(iim*3*jjm,llm) 
    193     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
    194     REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function 
    195  
    196     INTEGER :: i,j,ij,l 
    197     REAL(rstd) :: cp_ik, qv, temp, chi, nu, due, due_right, due_lup, due_ldown 
    198  
    199     CALL trace_start("compute_caldyn_fast") 
    200  
    201     IF(dysl_caldyn_fast) THEN 
    202 #include "../kernels_hex/caldyn_fast.k90" 
    203     ELSE 
    204  
    205     ! Compute Bernoulli term 
    206     IF(boussinesq) THEN 
    207        DO l=ll_begin,ll_end 
    208           !DIR$ SIMD 
    209           DO ij=ij_begin,ij_end          
    210              berni(ij,l) = pk(ij,l) 
    211              ! from now on pk contains the vertically-averaged geopotential 
    212              pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 
    213           END DO 
    214        END DO 
    215     ELSE ! compressible 
    216  
    217        DO l=ll_begin,ll_end 
    218           SELECT CASE(caldyn_thermo) 
    219           CASE(thermo_theta) ! vdp = theta.dpi => B = Phi 
    220              !DIR$ SIMD 
    221              DO ij=ij_begin,ij_end          
    222                 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 
    223              END DO 
    224           CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s)  
    225              !DIR$ SIMD 
    226              DO ij=ij_begin,ij_end          
    227                 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
    228                      + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy 
    229              END DO 
    230           CASE(thermo_moist)  
    231              !DIR$ SIMD 
    232              DO ij=ij_begin,ij_end          
    233                 ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T) 
    234                 ! Bd = Phi + gibbs_d 
    235                 ! Bv = Phi + gibbs_v 
    236                 ! pk=temperature, theta=entropy 
    237                 qv = theta(ij,l,2) 
    238                 temp = pk(ij,l) 
    239                 chi = log(temp/Treff) 
    240                 nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff) 
    241                 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
    242                      + temp*(cpp*(1.-chi)+Rd*nu) 
    243                 berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
    244                      + temp*(cppv*(1.-chi)+Rv*nu) 
    245             END DO 
    246           END SELECT 
    247        END DO 
    248  
    249     END IF ! Boussinesq/compressible 
    250  
    251 !!! u:=u+tau*du, du = -grad(B)-theta.grad(pi) 
    252     DO l=ll_begin,ll_end 
    253        IF(caldyn_thermo == thermo_moist) THEN 
    254           !DIR$ SIMD 
    255           DO ij=ij_begin,ij_end          
    256              due_right =  berni(ij+t_right,l)-berni(ij,l)       & 
    257                   + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   & 
    258                        *(pk(ij+t_right,l)-pk(ij,l))             & 
    259                   + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   & 
    260                        *(berniv(ij+t_right,l)-berniv(ij,l)) 
    261               
    262              due_lup = berni(ij+t_lup,l)-berni(ij,l)            & 
    263                   + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     & 
    264                        *(pk(ij+t_lup,l)-pk(ij,l))               & 
    265                   + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     & 
    266                        *(berniv(ij+t_lup,l)-berniv(ij,l)) 
    267               
    268              due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        & 
    269                   + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   & 
    270                        *(pk(ij+t_ldown,l)-pk(ij,l))             & 
    271                   + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   & 
    272                        *(berniv(ij+t_ldown,l)-berniv(ij,l)) 
    273               
    274              du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right 
    275              du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup 
    276              du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown 
    277              u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l) 
    278              u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l) 
    279              u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l) 
    280           END DO 
    281        ELSE 
    282           !DIR$ SIMD 
    283           DO ij=ij_begin,ij_end          
    284              due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  & 
    285                   *(pk(ij+t_right,l)-pk(ij,l))        & 
    286                   +  berni(ij+t_right,l)-berni(ij,l) 
    287              due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    & 
    288                   *(pk(ij+t_lup,l)-pk(ij,l))          & 
    289                   +  berni(ij+t_lup,l)-berni(ij,l) 
    290              due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) & 
    291                   *(pk(ij+t_ldown,l)-pk(ij,l))      & 
    292                   +   berni(ij+t_ldown,l)-berni(ij,l) 
    293              du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right 
    294              du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup 
    295              du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown 
    296              u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l) 
    297              u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l) 
    298              u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l) 
    299           END DO 
    300        END IF 
    301     END DO 
    302  
    303     END IF ! dysl 
    304     CALL trace_end("compute_caldyn_fast") 
    305  
    306   END SUBROUTINE compute_caldyn_fast 
    307  
    308185END MODULE caldyn_kernels_hevi_mod 
Note: See TracChangeset for help on using the changeset viewer.