Changeset 859


Ignore:
Timestamp:
05/07/19 01:09:27 (5 years ago)
Author:
jisesh
Message:

devel: separate module for compute_NH_geopot

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

Legend:

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

    r854 r859  
    1212  USE compute_caldyn_solver_mod, ONLY : compute_caldyn_solver 
    1313  USE compute_caldyn_fast_mod, ONLY : compute_caldyn_fast 
     14  USE compute_NH_geopot_mod, ONLY : compute_NH_geopot 
    1415  IMPLICIT NONE 
    1516  PRIVATE 
  • codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90

    r854 r859  
    99  PRIVATE 
    1010 
    11   REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6 
    12  
    13   LOGICAL, SAVE :: debug_hevi_solver = .FALSE. 
    14  
    15   PUBLIC :: compute_NH_geopot 
    1611 
    1712CONTAINS 
    1813 
    19   SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il) 
    20     REAL(rstd),INTENT(IN)    :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs 
    21     REAL(rstd),INTENT(IN)    :: phis(iim*jjm) 
    22     REAL(rstd),INTENT(IN)    :: m_ik(iim*jjm,llm) 
    23     REAL(rstd),INTENT(IN)    :: m_il(iim*jjm,llm+1) 
    24     REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm) 
    25     REAL(rstd),INTENT(IN)    :: W_il(iim*jjm,llm+1)   ! vertical momentum 
    26     REAL(rstd),INTENT(INOUT) :: Phi_il(iim*jjm,llm+1) ! geopotential 
    27  
    28     REAL(rstd) :: Phi_star_il(iim*jjm,llm+1)  
    29     REAL(rstd) :: p_ik(iim*jjm,llm)      ! pressure 
    30     REAL(rstd) :: R_il(iim*jjm,llm+1)    ! rhs of tridiag problem 
    31     REAL(rstd) :: x_il(iim*jjm,llm+1)    ! solution of tridiag problem 
    32     REAL(rstd) :: A_ik(iim*jjm,llm)      ! off-diagonal coefficients of tridiag problem 
    33     REAL(rstd) :: B_il(iim*jjm,llm+1)    ! diagonal coefficients of tridiag problem 
    34     REAL(rstd) :: C_ik(iim*jjm,llm)      ! Thomas algorithm 
    35     REAL(rstd) :: D_il(iim*jjm,llm+1)    ! Thomas algorithm 
    36     REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij 
    37     REAL(rstd) :: wil, tau2_g, g2, gm2, ml_g2, c2_mik, vreff 
    38  
    39     INTEGER    :: iter, ij, l, ij_omp_begin_ext, ij_omp_end_ext 
    40  
    41     CALL distrib_level(ij_begin_ext,ij_end_ext, ij_omp_begin_ext,ij_omp_end_ext) 
    42  
    43     IF(dysl) THEN 
    44 #define PHI_BOT(ij) phis(ij) 
    45 #include "../kernels_hex/compute_NH_geopot.k90" 
    46 #undef PHI_BOT 
    47     ELSE 
    48 !    FIXME : vertical OpenMP parallelism will not work 
    49     
    50     tau2_g=tau*tau/g 
    51     g2=g*g 
    52     gm2 = g**-2 
    53     gamma = 1./(1.-kappa) 
    54      
    55     ! compute Phi_star 
    56     DO l=1,llm+1 
    57        !DIR$ SIMD 
    58        DO ij=ij_begin_ext,ij_end_ext 
    59           Phi_star_il(ij,l) = Phi_il(ij,l) + tau*g2*(W_il(ij,l)/m_il(ij,l)-tau) 
    60        ENDDO 
    61     ENDDO 
    62  
    63     ! Newton-Raphson iteration : Phi_il contains current guess value 
    64     DO iter=1,5 ! 2 iterations should be enough 
    65  
    66        ! Compute pressure, A_ik 
    67        DO l=1,llm 
    68           !DIR$ SIMD 
    69           DO ij=ij_begin_ext,ij_end_ext 
    70              rho_ij = (g*m_ik(ij,l))/(Phi_il(ij,l+1)-Phi_il(ij,l)) 
    71              X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij 
    72              p_ik(ij,l) = preff*(X_ij**gamma) 
    73              c2_mik = gamma*p_ik(ij,l)/(rho_ij*m_ik(ij,l)) ! c^2 = gamma*R*T = gamma*p/rho 
    74              A_ik(ij,l) = c2_mik*(tau/g*rho_ij)**2 
    75           ENDDO 
    76        ENDDO 
    77  
    78        ! Compute residual, B_il 
    79        ! bottom interface l=1 
    80        !DIR$ SIMD 
    81        DO ij=ij_begin_ext,ij_end_ext 
    82           ml_g2 = gm2*m_il(ij,1)  
    83           B_il(ij,1) = A_ik(ij,1) + ml_g2 + tau2_g*rho_bot 
    84           R_il(ij,1) = ml_g2*( Phi_il(ij,1)-Phi_star_il(ij,1))  & 
    85                + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-phis(ij)) ) 
    86        ENDDO 
    87        ! inner interfaces 
    88        DO l=2,llm 
    89           !DIR$ SIMD 
    90           DO ij=ij_begin_ext,ij_end_ext 
    91              ml_g2 = gm2*m_il(ij,l)  
    92              B_il(ij,l) = A_ik(ij,l)+A_ik(ij,l-1) + ml_g2 
    93              R_il(ij,l) = ml_g2*( Phi_il(ij,l)-Phi_star_il(ij,l))  & 
    94                      + tau2_g*(p_ik(ij,l)-p_ik(ij,l-1)) 
    95              ! consistency check : if Wil=0 and initial state is in hydrostatic balance 
    96              ! then Phi_star_il(ij,l) = Phi_il(ij,l) - tau^2*g^2 
    97              ! and residual = tau^2*(ml+(1/g)dl_pi)=0 
    98           ENDDO 
    99        ENDDO 
    100        ! top interface l=llm+1 
    101        !DIR$ SIMD 
    102        DO ij=ij_begin_ext,ij_end_ext 
    103           ml_g2 = gm2*m_il(ij,llm+1)  
    104           B_il(ij,llm+1) = A_ik(ij,llm) + ml_g2 
    105           R_il(ij,llm+1) = ml_g2*( Phi_il(ij,llm+1)-Phi_star_il(ij,llm+1))  & 
    106                + tau2_g*( ptop-p_ik(ij,llm) ) 
    107        ENDDO 
    108  
    109        ! FIXME later 
    110        ! the lines below modify the tridiag problem  
    111        ! for flat, rigid boundary conditions at top and bottom : 
    112        ! zero out A(1), A(llm), R(1), R(llm+1) 
    113        ! => x(l)=0 at l=1,llm+1 
    114        DO ij=ij_begin_ext,ij_end_ext 
    115           A_ik(ij,1) = 0. 
    116           A_ik(ij,llm) = 0. 
    117           R_il(ij,1) = 0. 
    118           R_il(ij,llm+1) = 0. 
    119        ENDDO 
    120  
    121        IF(debug_hevi_solver) THEN ! print Linf(residual) 
    122           PRINT *, '[hevi_solver] R,p', iter, MAXVAL(ABS(R_il)), MAXVAL(p_ik) 
    123        END IF 
    124  
    125        ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm 
    126        ! Forward sweep : 
    127        ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),  
    128        ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1)) 
    129        ! bottom interface l=1 
    130        !DIR$ SIMD 
    131        DO ij=ij_begin_ext,ij_end_ext 
    132           X_ij = 1./B_il(ij,1) 
    133           C_ik(ij,1) = -A_ik(ij,1) * X_ij  
    134           D_il(ij,1) =  R_il(ij,1) * X_ij 
    135        ENDDO 
    136        ! inner interfaces/layers 
    137        DO l=2,llm 
    138           !DIR$ SIMD 
    139           DO ij=ij_begin_ext,ij_end_ext 
    140              X_ij = 1./(B_il(ij,l) + A_ik(ij,l-1)*C_ik(ij,l-1)) 
    141              C_ik(ij,l) = -A_ik(ij,l) * X_ij  
    142              D_il(ij,l) =  (R_il(ij,l)+A_ik(ij,l-1)*D_il(ij,l-1)) * X_ij 
    143           ENDDO 
    144        ENDDO 
    145        ! top interface l=llm+1 
    146        !DIR$ SIMD 
    147        DO ij=ij_begin_ext,ij_end_ext 
    148           X_ij = 1./(B_il(ij,llm+1) + A_ik(ij,llm)*C_ik(ij,llm)) 
    149           D_il(ij,llm+1) =  (R_il(ij,llm+1)+A_ik(ij,llm)*D_il(ij,llm)) * X_ij 
    150        ENDDO 
    151         
    152        ! Back substitution : 
    153        ! x(i) = D(i)-C(i)x(i+1), x(N+1)=0 
    154        ! + Newton-Raphson update 
    155        x_il=0. ! FIXME 
    156        ! top interface l=llm+1 
    157        !DIR$ SIMD 
    158        DO ij=ij_begin_ext,ij_end_ext 
    159           x_il(ij,llm+1) = D_il(ij,llm+1) 
    160           Phi_il(ij,llm+1) = Phi_il(ij,llm+1) - x_il(ij,llm+1) 
    161        ENDDO 
    162        ! lower interfaces 
    163        DO l=llm,1,-1 
    164           !DIR$ SIMD 
    165           DO ij=ij_begin_ext,ij_end_ext 
    166              x_il(ij,l) = D_il(ij,l) - C_ik(ij,l)*x_il(ij,l+1) 
    167              Phi_il(ij,l) = Phi_il(ij,l) - x_il(ij,l) 
    168           ENDDO 
    169        ENDDO 
    170  
    171        IF(debug_hevi_solver) THEN 
    172           PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il)) 
    173           PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il)) 
    174           DO l=1,llm+1 
    175              WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x', iter,l, MAXVAL(ABS(x_il(:,l))) 
    176           END DO 
    177        END IF 
    178  
    179     END DO ! Newton-Raphson 
    180  
    181     END IF ! dysl 
    182      
    183   END SUBROUTINE compute_NH_geopot 
    18414 
    18515END MODULE caldyn_kernels_hevi_mod 
  • codes/icosagcm/devel/src/dynamics/compute_caldyn_solver.F90

    r853 r859  
    1515    USE disvert_mod, ONLY : ptop 
    1616    USE caldyn_kernels_hevi_mod 
     17    USE compute_NH_geopot_mod 
    1718    REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6 
    1819    REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs 
Note: See TracChangeset for help on using the changeset viewer.