Changeset 878


Ignore:
Timestamp:
05/29/19 20:33:00 (5 years ago)
Author:
jisesh
Message:

devel: moved DYSL into compute_caldyn_solver.F90,compute_theta.F90 and compute_NH_geopot.F90

Location:
codes/icosagcm/devel
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/devel/Python/src/kernels_caldyn_NH.jin

    r876 r878  
    1010{%- endmacro %} 
    1111 
    12 KERNEL(compute_NH_geopot) 
    13   tau2_g=tau*tau/g 
    14   g2=g*g 
    15   gm2 = 1./g2 
    16   vreff = Treff*cpp/preff*kappa 
    17   gamma = 1./(1.-kappa) 
    18      
    19   BARRIER 
    20  
    21   ! compute Phi_star 
    22   SEQUENCE_C1 
    23     BODY('1,llm+1') 
    24       Phi_star_il(CELL) = Phi_il(CELL) + tau*g2*(W_il(CELL)/m_il(CELL)-tau) 
    25     END_BLOCK 
    26   END_BLOCK 
    27  
    28   ! Newton-Raphson iteration : Phi_il contains current guess value 
    29   DO iter=1,2 ! 2 iterations should be enough 
    30     ! Compute pressure, A_ik 
    31     SELECT CASE(caldyn_thermo) 
    32     CASE(thermo_theta) 
    33       {% call() compute_p_and_Aik() %} 
    34         X_ij = (cpp/preff)*kappa*theta(CELL)*rho_ij 
    35         p_ik(CELL) = preff*(X_ij**gamma) 
    36         c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho 
    37       {% endcall %} 
    38     CASE(thermo_entropy) 
    39       {% call() compute_p_and_Aik() %} 
    40         X_ij = log(vreff*rho_ij) + theta(CELL)/cpp 
    41         p_ik(CELL) = preff*exp(X_ij*gamma) 
    42         c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho 
    43       {% endcall %} 
    44     CASE DEFAULT 
    45         PRINT *, 'caldyn_thermo not supported by compute_NH_geopot', caldyn_thermo 
    46         STOP 
    47     END SELECT 
    48      
    49     ! NB : A(1), A(llm), R(1), R(llm+1) = 0 => x(l)=0 at l=1,llm+1 => flat, rigid top and bottom 
    50     ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm 
    51  
    52     SEQUENCE_C1 
    53       ! Compute residual R_il and B_il 
    54       PROLOGUE(1) 
    55         ! bottom interface l=1 
    56         ml_g2 = gm2*m_il(CELL)  
    57         B_il(CELL) = A_ik(CELL) + ml_g2 + tau2_g*rho_bot 
    58         R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  & 
    59                   + tau2_g*( p_ik(CELL)-pbot+rho_bot*(Phi_il(CELL)-PHI_BOT(HIDX(CELL))) ) 
    60       END_BLOCK 
    61       BODY('2,llm') 
    62         ! inner interfaces 
    63         ml_g2 = gm2*m_il(CELL)  
    64         B_il(CELL) = A_ik(CELL)+A_ik(DOWN(CELL)) + ml_g2 
    65         R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  & 
    66                   + tau2_g*(p_ik(CELL)-p_ik(DOWN(CELL))) 
    67         ! consistency check : if Wil=0 and initial state is in hydrostatic balance 
    68         ! then Phi_star_il(CELL) = Phi_il(CELL) - tau^2*g^2 
    69         ! and residual = tau^2*(ml+(1/g)dl_pi)=0 
    70       END_BLOCK 
    71       EPILOGUE('llm+1') 
    72         ! top interface l=llm+1 
    73         ml_g2 = gm2*m_il(CELL)  
    74         B_il(CELL) = A_ik(DOWN(CELL)) + ml_g2 
    75         R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  & 
    76                   + tau2_g*( ptop-p_ik(DOWN(CELL)) ) 
    77       END_BLOCK 
    78       ! 
    79       ! Forward sweep : 
    80       ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),  
    81       ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1)) 
    82       PROLOGUE(1) 
    83         X_ij = 1./B_il(CELL) 
    84         C_ik(CELL) = -A_ik(CELL) * X_ij  
    85         D_il(CELL) =  R_il(CELL) * X_ij 
    86       END_BLOCK 
    87       BODY('2,llm') 
    88         X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) ) 
    89         C_ik(CELL) = -A_ik(CELL) * X_ij  
    90         D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij 
    91       END_BLOCK 
    92       EPILOGUE('llm+1') 
    93         X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) ) 
    94         D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij 
    95         ! Back substitution : 
    96         ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=0 
    97         ! + Newton-Raphson update 
    98         ! top interface l=llm+1 
    99         x_il(CELL) = D_il(CELL) 
    100         Phi_il(CELL) = Phi_il(CELL) - x_il(CELL) 
    101       END_BLOCK 
    102       BODY('llm,1,-1') 
    103         ! Back substitution at lower interfaces 
    104         x_il(CELL) = D_il(CELL) - C_ik(CELL)*x_il(UP(CELL)) 
    105         Phi_il(CELL) = Phi_il(CELL) - x_il(CELL) 
    106       END_BLOCK 
    107     END_BLOCK 
    108  
    109     IF(debug_hevi_solver) THEN 
    110        PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il)) 
    111        PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il)) 
    112        DO l=1,llm+1 
    113           WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x_il', iter,l, MAXVAL(ABS(x_il(l,:))) 
    114        END DO 
    115        DO l=1,llm+1 
    116           WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] R_il', iter,l, MAXVAL(ABS(R_il(l,:)))  
    117        END DO 
    118     END IF 
    119  
    120   END DO ! Newton-Raphson 
    121  
    122   BARRIER 
    123  
    124   debug_hevi_solver=.FALSE. 
    125 END_BLOCK 
    126  
    12712KERNEL(caldyn_mil) 
    12813  FORALL_CELLS_EXT('1', 'llm+1') 
    … …  
    13116       CST_IF(IS_INNER_INTERFACE, m_il(CELL) = .5*(rhodz(CELL)+rhodz(DOWN(CELL))) ) 
    13217       CST_IF(IS_TOP_INTERFACE,   m_il(CELL) = .5*rhodz(DOWN(CELL))               ) 
    133     END_BLOCK 
    134   END_BLOCK 
    135 END_BLOCK 
    136  
    137 KERNEL(caldyn_solver) 
    138   !  
    139   ! Compute pressure (pres) and Exner function (pk) 
    140   ! kappa = R/Cp 
    141   ! 1-kappa = Cv/Cp 
    142   ! Cp/Cv = 1/(1-kappa) 
    143   gamma    = 1./(1.-kappa) 
    144   vreff    = Rd*Treff/preff ! reference specific volume 
    145   Cvd      = 1./(cpp-Rd) 
    146   Rd_preff = kappa*cpp/preff 
    147   FORALL_CELLS_EXT() 
    148     SELECT CASE(caldyn_thermo) 
    149     CASE(thermo_theta) 
    150       ON_PRIMAL 
    151         rho_ij = 1./(geopot(UP(CELL))-geopot(CELL)) 
    152         rho_ij = rho_ij*g*rhodz(CELL) 
    153         X_ij = Rd_preff*theta(CELL,1)*rho_ij 
    154         ! kappa.theta.rho = p/exner 
    155         ! => X = (p/p0)/(exner/Cp) 
    156         !      = (p/p0)^(1-kappa) 
    157         pres(CELL) = preff*(X_ij**gamma)  ! pressure 
    158         ! Compute Exner function (needed by compute_caldyn_fast) 
    159         ! other formulae possible if exponentiation is slow 
    160         pk(CELL)   = cpp*((pres(CELL)/preff)**kappa) ! Exner 
    161       END_BLOCK 
    162     CASE(thermo_entropy) 
    163       ON_PRIMAL 
    164         rho_ij = 1./(geopot(UP(CELL))-geopot(CELL)) 
    165         rho_ij = rho_ij*g*rhodz(CELL) 
    166         T_ij = Treff*exp( (theta(CELL,1)+Rd*log(vreff*rho_ij))*Cvd ) 
    167         pres(CELL) = rho_ij*Rd*T_ij 
    168         pk(CELL)   = T_ij 
    169       END_BLOCK 
    170     CASE DEFAULT 
    171       STOP 
    172     END SELECT 
    173   END_BLOCK 
    174  
    175   ! We need a barrier here because we compute pres above and do a vertical difference below 
    176   BARRIER 
    177  
    178   FORALL_CELLS_EXT('1', 'llm+1') 
    179     ON_PRIMAL 
    180       CST_IFTHEN(IS_BOTTOM_LEVEL) 
    181         ! Lower BC 
    182         dW(CELL)   = (1./g)*(pbot-rho_bot*(geopot(CELL)-PHI_BOT(HIDX(CELL)))-pres(CELL)) - m_il(CELL) 
    183       CST_ELSEIF(IS_TOP_INTERFACE) 
    184         ! Top BC 
    185         dW(CELL)   = (1./g)*(pres(DOWN(CELL))-ptop) - m_il(CELL) 
    186       CST_ELSE 
    187         dW(CELL)   = (1./g)*(pres(DOWN(CELL))-pres(CELL)) - m_il(CELL) 
    188       CST_ENDIF 
    189       W(CELL)    = W(CELL)+tau*dW(CELL) ! update W 
    190       dPhi(CELL) = g*g*W(CELL)/m_il(CELL) 
    191     END_BLOCK 
    192   END_BLOCK 
    193  
    194   ! We need a barrier here because we update W above and do a vertical average below 
    195   BARRIER 
    196  
    197   FORALL_CELLS_EXT() 
    198     ON_PRIMAL 
    199       ! compute du = -0.5*g^2.grad(w^2) 
    200       berni(CELL) = (-.25*g*g)*((W(CELL)/m_il(CELL))**2 + (W(UP(CELL))/m_il(UP(CELL)))**2 ) 
    201     END_BLOCK 
    202   END_BLOCK 
    203   FORALL_CELLS_EXT() 
    204     ON_EDGES 
    205       du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2)) 
    20618    END_BLOCK 
    20719  END_BLOCK 

  • codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin

    r876 r878  
    9696END_BLOCK 
    9797 
    98 KERNEL(theta) 
    99   IF(caldyn_eta==eta_mass) THEN ! Compute mass 
    100     ! FIXME : here mass_col is computed from rhodz 
    101     ! so that the DOFs are the same whatever caldyn_eta 
    102     ! in DYNAMICO mass_col is prognosed rather than rhodz 
    103     SEQUENCE_C1 
    104       PROLOGUE(0) 
    105         mass_col(HIDX(CELL))=0. 
    106       END_BLOCK 
    107       BODY('1,llm') 
    108           mass_col(HIDX(CELL)) = mass_col(HIDX(CELL)) + rhodz(CELL) 
    109       END_BLOCK 
    110     END_BLOCK 
    111     FORALL_CELLS_EXT() 
    112       ON_PRIMAL 
    113         ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on pressure rather than mass 
    114         !        m = mass_dak(CELL)+(mass_col(HIDX(CELL))*g+ptop)*mass_dbk(CELL) 
    115         !        rhodz(CELL) = m/g 
    116         rhodz(CELL) = mass_dak(CELL) + mass_col(HIDX(CELL))*mass_dbk(CELL) 
    117       END_BLOCK 
    118     END_BLOCK 
    119   END IF 
    120   DO iq=1,nqdyn 
    121     FORALL_CELLS_EXT() 
    122       ON_PRIMAL 
    123         theta(CELL,iq) = theta_rhodz(CELL,iq)/rhodz(CELL) 
    124       END_BLOCK 
    125     END_BLOCK 
    126   END DO 
    127 END_BLOCK 
    12898 

  • codes/icosagcm/devel/src/dynamics/compute_NH_geopot.F90

    r859 r878  
    66  LOGICAL, SAVE :: debug_hevi_solver = .FALSE. 
    77 
    8   PUBLIC :: compute_NH_geopot 
     8#include "../unstructured/unstructured.h90" 
     9 
     10  PUBLIC :: compute_NH_geopot,compute_NH_geopot_unst 
    911 
    1012CONTAINS 
     13 
     14#ifdef BEGIN_DYSL 
     15 
     16KERNEL(compute_NH_geopot) 
     17  tau2_g=tau*tau/g 
     18  g2=g*g 
     19  gm2 = 1./g2 
     20  vreff = Treff*cpp/preff*kappa 
     21  gamma = 1./(1.-kappa) 
     22 
     23  BARRIER 
     24 
     25  ! compute Phi_star 
     26  SEQUENCE_C1 
     27    BODY('1,llm+1') 
     28      Phi_star_il(CELL) = Phi_il(CELL) + tau*g2*(W_il(CELL)/m_il(CELL)-tau) 
     29    END_BLOCK 
     30  END_BLOCK 
     31 
     32  ! Newton-Raphson iteration : Phi_il contains current guess value 
     33  DO iter=1,2 ! 2 iterations should be enough 
     34    ! Compute pressure, A_ik 
     35    SELECT CASE(caldyn_thermo) 
     36    CASE(thermo_theta) 
     37      {% call() compute_p_and_Aik() %} 
     38        X_ij = (cpp/preff)*kappa*theta(CELL)*rho_ij 
     39        p_ik(CELL) = preff*(X_ij**gamma) 
     40        c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho 
     41      {% endcall %} 
     42    CASE(thermo_entropy) 
     43      {% call() compute_p_and_Aik() %} 
     44        X_ij = log(vreff*rho_ij) + theta(CELL)/cpp 
     45        p_ik(CELL) = preff*exp(X_ij*gamma) 
     46        c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho 
     47      {% endcall %} 
     48    CASE DEFAULT 
     49        PRINT *, 'caldyn_thermo not supported by compute_NH_geopot', caldyn_thermo 
     50        STOP 
     51    END SELECT 
     52 
     53    ! NB : A(1), A(llm), R(1), R(llm+1) = 0 => x(l)=0 at l=1,llm+1 => flat, rigid top and bottom 
     54    ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm 
     55 
     56    SEQUENCE_C1 
     57      ! Compute residual R_il and B_il 
     58      PROLOGUE(1) 
     59        ! bottom interface l=1 
     60        ml_g2 = gm2*m_il(CELL) 
     61        B_il(CELL) = A_ik(CELL) + ml_g2 + tau2_g*rho_bot 
     62        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  & 
     63                  + tau2_g*( p_ik(CELL)-pbot+rho_bot*(Phi_il(CELL)-PHI_BOT(HIDX(CELL))) ) 
     64      END_BLOCK 
     65      BODY('2,llm') 
     66        ! inner interfaces 
     67        ml_g2 = gm2*m_il(CELL) 
     68        B_il(CELL) = A_ik(CELL)+A_ik(DOWN(CELL)) + ml_g2 
     69        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  & 
     70                  + tau2_g*(p_ik(CELL)-p_ik(DOWN(CELL))) 
     71        ! consistency check : if Wil=0 and initial state is in hydrostatic balance 
     72        ! then Phi_star_il(CELL) = Phi_il(CELL) - tau^2*g^2 
     73        ! and residual = tau^2*(ml+(1/g)dl_pi)=0 
     74      END_BLOCK 
     75      EPILOGUE('llm+1') 
     76        ! top interface l=llm+1 
     77        ml_g2 = gm2*m_il(CELL) 
     78        B_il(CELL) = A_ik(DOWN(CELL)) + ml_g2 
     79        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  & 
     80                  + tau2_g*( ptop-p_ik(DOWN(CELL)) ) 
     81      END_BLOCK 
     82      ! 
     83      ! Forward sweep : 
     84      ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)), 
     85      ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1)) 
     86      PROLOGUE(1) 
     87        X_ij = 1./B_il(CELL) 
     88        C_ik(CELL) = -A_ik(CELL) * X_ij 
     89        D_il(CELL) =  R_il(CELL) * X_ij 
     90      END_BLOCK 
     91      BODY('2,llm') 
     92        X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) ) 
     93        C_ik(CELL) = -A_ik(CELL) * X_ij 
     94        D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij 
     95      END_BLOCK 
     96      EPILOGUE('llm+1') 
     97        X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) ) 
     98        D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij 
     99        ! Back substitution : 
     100        ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=0 
     101        ! + Newton-Raphson update 
     102        ! top interface l=llm+1 
     103        x_il(CELL) = D_il(CELL) 
     104        Phi_il(CELL) = Phi_il(CELL) - x_il(CELL) 
     105      END_BLOCK 
     106      BODY('llm,1,-1') 
     107        ! Back substitution at lower interfaces 
     108        x_il(CELL) = D_il(CELL) - C_ik(CELL)*x_il(UP(CELL)) 
     109        Phi_il(CELL) = Phi_il(CELL) - x_il(CELL) 
     110      END_BLOCK 
     111    END_BLOCK 
     112 
     113    IF(debug_hevi_solver) THEN 
     114       PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il)) 
     115       PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il)) 
     116       DO l=1,llm+1 
     117          WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x_il', iter,l, MAXVAL(ABS(x_il(l,:))) 
     118       END DO 
     119       DO l=1,llm+1 
     120          WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] R_il', iter,l, MAXVAL(ABS(R_il(l,:))) 
     121       END DO 
     122    END IF 
     123 
     124  END DO ! Newton-Raphson 
     125 
     126  BARRIER 
     127 
     128  debug_hevi_solver=.FALSE. 
     129END_BLOCK 
     130 
     131#endif END_DYSL 
     132 
     133SUBROUTINE compute_NH_geopot_unst(tau, m_ik, m_il, theta, W_il, Phi_il) 
     134  USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT 
     135  USE disvert_mod, only : g,Treff,cpp,preff,kappa,caldyn_thermo,thermo_theta, & 
     136    thermo_entropy 
     137  USE disvert_mod, ONLY : ptop 
     138  USE data_unstructured_mod, ONLY : primal_num,edge_num,dual_num,id_NH_geopot,debug_hevi_solver_, & 
     139    PHI_BOT,pbot,rho_bot,enter_trace, exit_trace 
     140  FIELD_MASS   :: m_ik, theta   ! IN*2 
     141  FIELD_GEOPOT :: m_il, W_il, Phi_il, Phi_star_il  ! IN,INOUT*2, LOCAL 
     142  NUM :: tau, gamma, tau2_g, tau2_g2, g2, gm2, vreff, Rd_preff 
     143  INTEGER :: iter 
     144  LOGICAL :: debug_hevi_solver 
     145  DECLARE_INDICES 
     146  NUM :: rho_ij, X_ij, Y_ij, wil, rho_c2_mik, c2_mik, ml_g2 
     147#define COLUMN 0 
     148#if COLUMN  
     149  NUM1(llm)  :: pk, Ak, Ck 
     150  NUM1(llm+1):: Rl, Bl, Dl, xl 
     151#define p_ik(l,ij) pk(l) 
     152#define A_ik(l,ij) Ak(l) 
     153#define C_ik(l,ij) Ck(l) 
     154#define R_il(l,ij) Rl(l) 
     155#define B_il(l,ij) Bl(l) 
     156#define D_il(l,ij) Dl(l) 
     157#define x_il(l,ij) xl(l) 
     158#else 
     159  FIELD_MASS :: p_ik, A_ik, C_ik 
     160  FIELD_GEOPOT :: R_il, B_il, D_il, x_il 
     161#endif 
     162 
     163  debug_hevi_solver=.FALSE. 
     164!$OMP MASTER 
     165  debug_hevi_solver = debug_hevi_solver_ 
     166!$OMP END MASTER 
     167 
     168#define PHI_BOT(ij) Phi_bot 
     169  START_TRACE(id_NH_geopot, 7,0,0) 
     170#include "../kernels_unst/compute_NH_geopot.k90" 
     171  STOP_TRACE 
     172 
     173!$OMP MASTER 
     174  debug_hevi_solver_ = debug_hevi_solver 
     175!$OMP END MASTER 
     176#undef PHI_BOT 
     177 
     178#if COLUMN 
     179#undef p_ik 
     180#undef A_ik 
     181#undef C_ik 
     182#undef R_il 
     183#undef B_il 
     184#undef D_il 
     185#undef x_il 
     186#endif 
     187#undef COLUMN 
     188END SUBROUTINE compute_NH_geopot_unst 
    11189 
    12190  SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il) 

  • codes/icosagcm/devel/src/dynamics/compute_caldyn_solver.F90

    r859 r878  
    44  PRIVATE 
    55 
     6#include "../unstructured/unstructured.h90" 
     7 
    68  PUBLIC :: compute_caldyn_solver 
    79 
    810CONTAINS 
     11 
     12#ifdef BEGIN_DYSL 
     13 
     14KERNEL(caldyn_solver) 
     15  !    
     16  ! Compute pressure (pres) and Exner function (pk) 
     17  ! kappa = R/Cp 
     18  ! 1-kappa = Cv/Cp 
     19  ! Cp/Cv = 1/(1-kappa) 
     20  gamma    = 1./(1.-kappa) 
     21  vreff    = Rd*Treff/preff ! reference specific volume 
     22  Cvd      = 1./(cpp-Rd) 
     23  Rd_preff = kappa*cpp/preff 
     24  FORALL_CELLS_EXT() 
     25    SELECT CASE(caldyn_thermo) 
     26    CASE(thermo_theta) 
     27      ON_PRIMAL 
     28        rho_ij = 1./(geopot(UP(CELL))-geopot(CELL)) 
     29        rho_ij = rho_ij*g*rhodz(CELL) 
     30        X_ij = Rd_preff*theta(CELL,1)*rho_ij 
     31        ! kappa.theta.rho = p/exner 
     32        ! => X = (p/p0)/(exner/Cp) 
     33        !      = (p/p0)^(1-kappa) 
     34        pres(CELL) = preff*(X_ij**gamma)  ! pressure 
     35        ! Compute Exner function (needed by compute_caldyn_fast) 
     36        ! other formulae possible if exponentiation is slow 
     37        pk(CELL)   = cpp*((pres(CELL)/preff)**kappa) ! Exner 
     38      END_BLOCK 
     39    CASE(thermo_entropy) 
     40      ON_PRIMAL 
     41        rho_ij = 1./(geopot(UP(CELL))-geopot(CELL)) 
     42        rho_ij = rho_ij*g*rhodz(CELL) 
     43        T_ij = Treff*exp( (theta(CELL,1)+Rd*log(vreff*rho_ij))*Cvd ) 
     44        pres(CELL) = rho_ij*Rd*T_ij 
     45        pk(CELL)   = T_ij 
     46      END_BLOCK 
     47    CASE DEFAULT 
     48      STOP 
     49    END SELECT 
     50  END_BLOCK 
     51 
     52  ! We need a barrier here because we compute pres above and do a vertical difference below 
     53  BARRIER 
     54 
     55  FORALL_CELLS_EXT('1', 'llm+1') 
     56    ON_PRIMAL 
     57      CST_IFTHEN(IS_BOTTOM_LEVEL) 
     58        ! Lower BC 
     59        dW(CELL)   = (1./g)*(pbot-rho_bot*(geopot(CELL)-PHI_BOT(HIDX(CELL)))-pres(CELL)) - m_il(CELL) 
     60      CST_ELSEIF(IS_TOP_INTERFACE) 
     61        ! Top BC 
     62        dW(CELL)   = (1./g)*(pres(DOWN(CELL))-ptop) - m_il(CELL) 
     63      CST_ELSE 
     64        dW(CELL)   = (1./g)*(pres(DOWN(CELL))-pres(CELL)) - m_il(CELL) 
     65      CST_ENDIF 
     66      W(CELL)    = W(CELL)+tau*dW(CELL) ! update W 
     67      dPhi(CELL) = g*g*W(CELL)/m_il(CELL) 
     68    END_BLOCK 
     69  END_BLOCK 
     70 
     71  ! We need a barrier here because we update W above and do a vertical average below 
     72  BARRIER 
     73 
     74  FORALL_CELLS_EXT() 
     75    ON_PRIMAL 
     76      ! compute du = -0.5*g^2.grad(w^2) 
     77      berni(CELL) = (-.25*g*g)*((W(CELL)/m_il(CELL))**2 + (W(UP(CELL))/m_il(UP(CELL)))**2 ) 
     78    END_BLOCK 
     79  END_BLOCK 
     80  FORALL_CELLS_EXT() 
     81    ON_EDGES 
     82      du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2)) 
     83    END_BLOCK 
     84  END_BLOCK 
     85END_BLOCK 
     86 
     87#endif END_DYSL 
     88 
     89SUBROUTINE compute_caldyn_solver_unst(tau,rhodz,theta, berni,pres,m_il, pk,geopot,W,dPhi,dW,du) 
     90  USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT 
     91  USE earth_const 
     92  USE trace 
     93  USE grid_param, ONLY : nqdyn 
     94  USE disvert_mod, ONLY : ptop 
     95  USE data_unstructured_mod, ONLY : id_solver,primal_num,dual_num,edge_num,left, right,PHI_BOT, & 
     96    enter_trace, exit_trace 
     97  USE compute_NH_geopot_mod, ONLY : compute_NH_geopot_unst 
     98  NUM, INTENT(IN) :: tau  
     99  FIELD_MASS   :: rhodz,pk,berni,pres    ! IN, OUT, BUF*2 
     100  FIELD_THETA  :: theta                  ! IN 
     101  FIELD_GEOPOT :: geopot,W,dPhi,dW, m_il ! INOUT,INOUT, OUT,OUT, BUF 
     102  FIELD_U      :: du                     ! OUT 
     103  DECLARE_INDICES 
     104  NUM :: X_ij, rho_ij, T_ij, gamma, Cvd, vreff, Rd_preff 
     105  REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6  ! FIXME 
     106#define PHI_BOT(ij) Phi_bot 
     107#include "../kernels_unst/caldyn_mil.k90" 
     108  IF(tau>0) THEN ! solve implicit problem for geopotential 
     109    CALL compute_NH_geopot_unst(tau, rhodz, m_il, theta, W, geopot) 
     110  END IF 
     111  START_TRACE(id_solver, 7,0,1) 
     112#include "../kernels_unst/caldyn_solver.k90" 
     113  STOP_TRACE 
     114#undef PHI_BOT 
     115END SUBROUTINE compute_caldyn_solver_unst 
    9116 
    10117  SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du) 

  • codes/icosagcm/devel/src/dynamics/compute_theta.F90

    r831 r878  
    44  PRIVATE 
    55 
     6#include "../unstructured/unstructured.h90" 
     7 
    68  PUBLIC :: compute_theta 
    79 
    810CONTAINS 
     11 
     12#ifdef BEGIN_DYSL 
     13 
     14KERNEL(theta) 
     15  IF(caldyn_eta==eta_mass) THEN ! Compute mass 
     16    ! FIXME : here mass_col is computed from rhodz 
     17    ! so that the DOFs are the same whatever caldyn_eta 
     18    ! in DYNAMICO mass_col is prognosed rather than rhodz 
     19    SEQUENCE_C1 
     20      PROLOGUE(0) 
     21        mass_col(HIDX(CELL))=0. 
     22      END_BLOCK 
     23      BODY('1,llm') 
     24          mass_col(HIDX(CELL)) = mass_col(HIDX(CELL)) + rhodz(CELL) 
     25      END_BLOCK 
     26    END_BLOCK 
     27    FORALL_CELLS_EXT() 
     28      ON_PRIMAL 
     29        ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on 
     30        ! pressure rather than mass 
     31        !        m = mass_dak(CELL)+(mass_col(HIDX(CELL))*g+ptop)*mass_dbk(CELL) 
     32        !        rhodz(CELL) = m/g  
     33        rhodz(CELL) = MASS_DAK(CELL) + mass_col(HIDX(CELL))*MASS_DBK(CELL) 
     34      END_BLOCK 
     35    END_BLOCK 
     36  END IF 
     37  DO iq=1,nqdyn 
     38    FORALL_CELLS_EXT() 
     39      ON_PRIMAL 
     40        theta(CELL,iq) = theta_rhodz(CELL,iq)/rhodz(CELL) 
     41      END_BLOCK 
     42    END_BLOCK 
     43  END DO 
     44END_BLOCK 
     45 
     46#endif END_DYSL 
     47 
     48SUBROUTINE compute_theta_unst(mass_col,rhodz,theta_rhodz, theta) 
     49  USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT 
     50  USE data_unstructured_mod, ONLY : id_theta,primal_num,dual_num,edge_num, & 
     51    enter_trace, exit_trace 
     52  USE grid_param, ONLY : nqdyn 
     53  USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop 
     54  FIELD_PS :: mass_col 
     55  FIELD_MASS :: rhodz 
     56  FIELD_THETA :: theta, theta_rhodz 
     57  DECLARE_INDICES 
     58  NUM :: m 
     59  START_TRACE(id_theta, 3,0,0) ! primal, dual, edge   
     60#define MASS_DAK(l,ij) mass_dak(l) 
     61#define MASS_DBK(l,ij) mass_dbk(l) 
     62#include "../kernels_unst/theta.k90" 
     63#undef MASS_DAK 
     64#undef MASS_DBK 
     65  STOP_TRACE 
     66END SUBROUTINE 
    967 
    1068  SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta) 

  • codes/icosagcm/devel/src/kernels_unst/theta.k90

    r686 r878  
    2020            ! m = mass_dak(l,ij)+(mass_col(ij)*g+ptop)*mass_dbk(l,ij) 
    2121            ! rhodz(l,ij) = m/g 
    22             rhodz(l,ij) = mass_dak(l,ij) + mass_col(ij)*mass_dbk(l,ij) 
     22            rhodz(l,ij) = MASS_DAK(l,ij) + mass_col(ij)*MASS_DBK(l,ij) 
    2323         END DO 
    2424      END DO 
Note: See TracChangeset for help on using the changeset viewer.