Changeset 368 for codes


Ignore:
Timestamp:
01/13/16 17:36:50 (8 years ago)
Author:
dubos
Message:

NH HEVI implicit solver - tested with DCMIP3

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r366 r368  
    66  PRIVATE 
    77 
     8  REAL(rstd), PARAMETER :: pbot=1e5, Phi_bot=0., rho_bot=1e6 ! FIXME 
     9 
     10  LOGICAL, PARAMETER :: debug_hevi_solver = .FALSE. 
     11 
    812  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_slow, & 
    913       compute_caldyn_solver, compute_caldyn_fast 
    … …  
    1216 
    1317  SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta) 
    14     USE icosa 
    1518    USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop 
    1619    USE exner_mod 
    … …  
    5255 
    5356  SUBROUTINE compute_pvort_only(u,rhodz,qu,qv) 
    54     USE icosa 
    5557    USE exner_mod 
    5658    USE trace 
    … …  
    121123  END SUBROUTINE compute_pvort_only 
    122124 
     125  SUBROUTINE compute_NH_geopot(tau, m_ik, m_il, theta, W_il, Phi_il) 
     126    USE omp_para 
     127    USE disvert_mod 
     128    IMPLICIT NONE 
     129    REAL(rstd),INTENT(IN)    :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs 
     130    REAL(rstd),INTENT(IN)    :: m_ik(iim*jjm,llm) 
     131    REAL(rstd),INTENT(IN)    :: m_il(iim*jjm,llm+1) 
     132    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm) 
     133    REAL(rstd),INTENT(IN)    :: W_il(iim*jjm,llm+1)   ! vertical momentum 
     134    REAL(rstd),INTENT(INOUT) :: Phi_il(iim*jjm,llm+1) ! geopotential 
     135 
     136    REAL(rstd) :: Phi_star_il(iim*jjm,llm+1)  
     137    REAL(rstd) :: p_ik(iim*jjm,llm)      ! pressure 
     138    REAL(rstd) :: R_il(iim*jjm,llm+1)    ! rhs of tridiag problem 
     139    REAL(rstd) :: x_il(iim*jjm,llm+1)    ! solution of tridiag problem 
     140    REAL(rstd) :: A_ik(iim*jjm,llm)      ! off-diagonal coefficients of tridiag problem 
     141    REAL(rstd) :: B_il(iim*jjm,llm+1)    ! diagonal coefficients of tridiag problem 
     142    REAL(rstd) :: C_ik(iim*jjm,llm)      ! Thomas algorithm 
     143    REAL(rstd) :: D_il(iim*jjm,llm+1)    ! Thomas algorithm 
     144    REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij 
     145    REAL(rstd) :: wil, tau2_g, g2, gm2, ml_g2, c2_mik 
     146 
     147    INTEGER    :: iter, ij, l 
     148 
     149!    FIXME : vertical OpenMP parallelism will not work 
     150    
     151    tau2_g=tau*tau/g 
     152    g2=g*g 
     153    gm2 = g**-2 
     154    gamma = 1./(1.-kappa) 
     155     
     156    ! compute Phi_star 
     157    DO l=1,llm+1 
     158       !DIR$ SIMD 
     159       DO ij=ij_begin_ext,ij_end_ext 
     160          Phi_star_il(ij,l) = Phi_il(ij,l) + tau*g2*(W_il(ij,l)/m_il(ij,l)-tau) 
     161       ENDDO 
     162    ENDDO 
     163 
     164    ! Newton-Raphson iteration : Phi_il contains current guess value 
     165    DO iter=1,2 ! 2 iterations should be enough 
     166 
     167       ! Compute pressure, A_ik 
     168       DO l=1,llm 
     169          !DIR$ SIMD 
     170          DO ij=ij_begin_ext,ij_end_ext 
     171             rho_ij = (g*m_ik(ij,l))/(Phi_il(ij,l+1)-Phi_il(ij,l)) 
     172             X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij 
     173             p_ik(ij,l) = preff*(X_ij**gamma) 
     174             c2_mik = gamma*p_ik(ij,l)/(rho_ij*m_ik(ij,l)) ! c^2 = gamma*R*T = gamma*p/rho 
     175             A_ik(ij,l) = c2_mik*(tau/g*rho_ij)**2 
     176          ENDDO 
     177       ENDDO 
     178 
     179       ! Compute residual, B_il 
     180       ! bottom interface l=1 
     181       !DIR$ SIMD 
     182       DO ij=ij_begin_ext,ij_end_ext 
     183          ml_g2 = gm2*m_il(ij,1)  
     184          B_il(ij,1) = A_ik(ij,1) + ml_g2 + tau2_g*rho_bot 
     185          R_il(ij,1) = ml_g2*( Phi_il(ij,1)-Phi_star_il(ij,1))  & 
     186               + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-Phi_bot) ) 
     187       ENDDO 
     188       ! inner interfaces 
     189       DO l=2,llm 
     190          !DIR$ SIMD 
     191          DO ij=ij_begin_ext,ij_end_ext 
     192             ml_g2 = gm2*m_il(ij,l)  
     193             B_il(ij,l) = A_ik(ij,l)+A_ik(ij,l-1) + ml_g2 
     194             R_il(ij,l) = ml_g2*( Phi_il(ij,l)-Phi_star_il(ij,l))  & 
     195                     + tau2_g*(p_ik(ij,l)-p_ik(ij,l-1)) 
     196             ! consistency check : if Wil=0 and initial state is in hydrostatic balance 
     197             ! then Phi_star_il(ij,l) = Phi_il(ij,l) - tau^2*g^2 
     198             ! and residual = tau^2*(ml+(1/g)dl_pi)=0 
     199          ENDDO 
     200       ENDDO 
     201       ! top interface l=llm+1 
     202       !DIR$ SIMD 
     203       DO ij=ij_begin_ext,ij_end_ext 
     204          ml_g2 = gm2*m_il(ij,llm+1)  
     205          B_il(ij,llm+1) = A_ik(ij,llm) + ml_g2 
     206          R_il(ij,llm+1) = ml_g2*( Phi_il(ij,llm+1)-Phi_star_il(ij,llm+1))  & 
     207               + tau2_g*( ptop-p_ik(ij,llm) ) 
     208       ENDDO 
     209 
     210       ! FIXME later 
     211       ! the lines below modify the tridiag problem  
     212       ! for flat, rigid boundary conditions at top and bottom : 
     213       ! zero out A(1), A(llm), R(1), R(llm+1) 
     214       ! => x(l)=0 at l=1,llm+1 
     215       DO ij=ij_begin_ext,ij_end_ext 
     216          A_ik(ij,1) = 0. 
     217          A_ik(ij,llm) = 0. 
     218          R_il(ij,1) = 0. 
     219          R_il(ij,llm+1) = 0. 
     220       ENDDO 
     221 
     222       IF(debug_hevi_solver) THEN ! print Linf(residual) 
     223          PRINT *, '[hevi_solver] R,p', iter, MAXVAL(ABS(R_il)), MAXVAL(p_ik) 
     224       END IF 
     225 
     226       ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm 
     227       ! Forward sweep : 
     228       ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),  
     229       ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1)) 
     230       ! bottom interface l=1 
     231       !DIR$ SIMD 
     232       DO ij=ij_begin_ext,ij_end_ext 
     233          X_ij = 1./B_il(ij,1) 
     234          C_ik(ij,1) = -A_ik(ij,1) * X_ij  
     235          D_il(ij,1) =  R_il(ij,1) * X_ij 
     236       ENDDO 
     237       ! inner interfaces/layers 
     238       DO l=2,llm 
     239          !DIR$ SIMD 
     240          DO ij=ij_begin_ext,ij_end_ext 
     241             X_ij = 1./(B_il(ij,l) + A_ik(ij,l-1)*C_ik(ij,l-1)) 
     242             C_ik(ij,l) = -A_ik(ij,l) * X_ij  
     243             D_il(ij,l) =  (R_il(ij,l)+A_ik(ij,l-1)*D_il(ij,l-1)) * X_ij 
     244          ENDDO 
     245       ENDDO 
     246       ! top interface l=llm+1 
     247       !DIR$ SIMD 
     248       DO ij=ij_begin_ext,ij_end_ext 
     249          X_ij = 1./(B_il(ij,llm+1) + A_ik(ij,llm)*C_ik(ij,llm)) 
     250          D_il(ij,llm+1) =  (R_il(ij,llm+1)+A_ik(ij,llm)*D_il(ij,llm)) * X_ij 
     251       ENDDO 
     252        
     253       ! Back substitution : 
     254       ! x(i) = D(i)-C(i)x(i+1), x(N+1)=0 
     255       ! + Newton-Raphson update 
     256       x_il=0. ! FIXME 
     257       ! top interface l=llm+1 
     258       !DIR$ SIMD 
     259       DO ij=ij_begin_ext,ij_end_ext 
     260          x_il(ij,llm+1) = D_il(ij,llm+1) 
     261          Phi_il(ij,llm+1) = Phi_il(ij,llm+1) - x_il(ij,llm+1) 
     262       ENDDO 
     263       ! lower interfaces 
     264       DO l=llm,1,-1 
     265          !DIR$ SIMD 
     266          DO ij=ij_begin_ext,ij_end_ext 
     267             x_il(ij,l) = D_il(ij,l) - C_ik(ij,l)*x_il(ij,l+1) 
     268             Phi_il(ij,l) = Phi_il(ij,l) - x_il(ij,l) 
     269          ENDDO 
     270       ENDDO 
     271 
     272       IF(debug_hevi_solver) THEN 
     273          PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il)) 
     274          PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il)) 
     275          DO l=1,llm+1 
     276             WRITE(*,'(A,I2.1,I3.2,E9.2)'), '[hevi_solver] x', iter,l, MAXVAL(ABS(x_il(:,l))) 
     277          END DO 
     278       END IF 
     279 
     280    END DO ! Newton-Raphson 
     281     
     282  END SUBROUTINE compute_NH_geopot 
     283 
    123284  SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk, geopot,W, dPhi,dW) 
    124     USE icosa 
    125285    USE disvert_mod 
    126286    USE exner_mod 
    127287    USE trace 
    128288    USE omp_para 
    129     USE disvert_mod 
    130289    IMPLICIT NONE 
    131290    REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs 
    … …  
    138297    REAL(rstd),INTENT(OUT)   :: dPhi(iim*jjm,llm+1) 
    139298 
    140     INTEGER :: ij,l 
    141     INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext 
    142     REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij, m_il 
     299    REAL(rstd) :: m_il(iim*jjm,llm+1)        ! rhodz averaged to interfaces 
     300    REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij 
     301    INTEGER    :: ij, l 
    143302 
    144303    CALL trace_start("compute_caldyn_solver") 
    145304 
    146     CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext) 
    147     ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1 
    148     ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1 
    149  
    150     IF(tau>0) THEN 
    151        PRINT *, 'No implicit solver for NH yet' 
    152        STOP 
     305    ! FIXME : vertical OpenMP parallelism will not work 
     306 
     307    ! average m_ik to interfaces => m_il 
     308    !DIR$ SIMD 
     309    DO ij=ij_begin_ext,ij_end_ext 
     310       m_il(ij,1) = .5*rhodz(ij,1) 
     311    ENDDO 
     312    DO l=2,llm 
     313       !DIR$ SIMD 
     314       DO ij=ij_begin_ext,ij_end_ext 
     315          m_il(ij,l) = .5*(rhodz(ij,l-1)+rhodz(ij,l)) 
     316       ENDDO 
     317    ENDDO 
     318    !DIR$ SIMD 
     319    DO ij=ij_begin_ext,ij_end_ext 
     320       m_il(ij,llm+1) = .5*rhodz(ij,llm) 
     321    ENDDO 
     322 
     323    IF(tau>0) THEN ! solve implicit problem for geopotential 
     324       CALL compute_NH_geopot(tau, rhodz, m_il, theta, W, geopot) 
    153325    END IF 
    154326 
    … …  
    158330    ! Cp/Cv = 1/(1-kappa) 
    159331    gamma = 1./(1.-kappa) 
    160     DO l=ll_begin,ll_end 
    161        !DIR$ SIMD 
    162        DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     332    DO l=1,llm 
     333       !DIR$ SIMD 
     334       DO ij=ij_begin_ext,ij_end_ext 
    163335          rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l)) 
    164336          X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij 
    … …  
    170342    ENDDO 
    171343 
    172     ! Compute tendencies for geopotential and vertical momentum 
    173     DO l=ll_beginp1,ll_end 
    174        !DIR$ SIMD 
    175        DO ij=ij_omp_begin_ext,ij_omp_end_ext 
    176           m_il = .5*(rhodz(ij,l)+rhodz(ij,l-1)) 
    177           dW(ij,l) = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il 
    178           dPhi(ij,l) = g*g*W(ij,l)/m_il 
     344    ! Update W, compute tendencies for geopotential and vertical momentum 
     345    DO l=2,llm 
     346       !DIR$ SIMD 
     347       DO ij=ij_begin_ext,ij_end_ext 
     348          dW(ij,l)   = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il(ij,l) 
     349          W(ij,l)    = W(ij,l)+tau*dW(ij,l) ! update W 
     350          dPhi(ij,l) = g*g*W(ij,l)/m_il(ij,l) 
    179351       ENDDO 
    180352       !       PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l))) 
    … …  
    191363    ENDDO 
    192364    ! Upper BC p=ptop 
    193 !    DO ij=ij_omp_begin_ext,ij_omp_end_ext 
    194 !       dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm) 
    195 !       dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm) 
    196 !    ENDDO 
     365    !    DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     366    !       dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm) 
     367    !       dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm) 
     368    !    ENDDO 
    197369     
    198370    ! Compute Exner function (needed by compute_caldyn_fast) 
    199     DO l=ll_begin,ll_end 
    200        !DIR$ SIMD 
    201        DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     371    DO l=1,llm 
     372       !DIR$ SIMD 
     373       DO ij=ij_begin_ext,ij_end_ext 
    202374          pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow 
    203375       ENDDO 
Note: See TracChangeset for help on using the changeset viewer.