Changeset 580


Ignore:
Timestamp:
10/13/17 15:55:07 (7 years ago)
Author:
dubos
Message:

trunk : upgrading to devel

Location:
codes/icosagcm/trunk/src
Files:
12 added
6 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/base/earth_const.f90

    r553 r580  
    2828  !$OMP THREADPRIVATE(boussinesq) 
    2929  !$OMP THREADPRIVATE(hydrostatic)  
     30  LOGICAL :: dysl, dysl_geopot, dysl_pvort_only, dysl_caldyn_fast, dysl_caldyn_coriolis, dysl_slow_hydro 
     31  !$OMP THREADPRIVATE(dysl, dysl_geopot, dysl_pvort_only, dysl_caldyn_fast, dysl_caldyn_coriolis, dysl_slow_hydro) 
    3032 
    3133CONTAINS 

  • codes/icosagcm/trunk/src/dynamics/caldyn.f90

    r548 r580  
    22  USE icosa 
    33  PRIVATE 
     4  SAVE 
    45  CHARACTER(LEN=255),SAVE :: caldyn_type 
    56!$OMP THREADPRIVATE(caldyn_type) 
    6    
    7   PUBLIC init_caldyn, caldyn, caldyn_BC 
     7  
     8  PUBLIC init_caldyn, caldyn, caldyn_BC, dysl 
    89   
    910CONTAINS 
    … …  
    2829    END SELECT 
    2930         
    30        
    3131  END SUBROUTINE init_caldyn 
    3232 

  • codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90

    r561 r580  
    2222    REAL(rstd),POINTER :: planetvel(:) 
    2323 
    24 #ifdef CPP_DYSL 
    25        IF(is_master) PRINT *,'CPP_DYSL : Using macro-generated compute kernels' 
    26 #endif 
    27  
    2824    hydrostatic=.TRUE. 
    2925    CALL getin("hydrostatic",hydrostatic) 
    3026   
     27    dysl=.NOT.hydrostatic ! dysl defaults to .FALSE. if hydrostatic, .TRUE. if NH                                               
     28    CALL getin("dysl",dysl) 
     29    dysl_geopot=dysl 
     30    CALL getin("dysl_geopot",dysl_geopot) 
     31    dysl_caldyn_fast=dysl 
     32    CALL getin("dysl_caldyn_fast",dysl_caldyn_fast) 
     33    dysl_slow_hydro=dysl 
     34    CALL getin("dysl_slow_hydro",dysl_slow_hydro) 
     35    dysl_pvort_only=dysl 
     36    CALL getin("dysl_pvort_only",dysl_pvort_only) 
     37    dysl_caldyn_coriolis=dysl 
     38    CALL getin("dysl_caldyn_coriolis",dysl_caldyn_coriolis) 
     39 
    3140    def='energy' 
    3241    CALL getin('caldyn_conserv',def) 

  • codes/icosagcm/trunk/src/dynamics/caldyn_hevi.f90

    r561 r580  
    4848    TYPE(t_field) :: f_dPhi_fast(:) 
    4949     
    50     REAL(rstd),POINTER :: ps(:), dps(:) 
     50    REAL(rstd),POINTER :: ps(:), dps(:), phis(:) 
    5151    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:) 
    5252    REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:) 
    … …  
    115115          CALL compute_geopot(mass,theta, ps,pk,geopot) 
    116116       ELSE 
     117          phis = f_phis(ind) 
    117118          W = f_W(ind) 
    118119          dW = f_dW_fast(ind) 
    … …  
    121122          m_il = f_wil(ind) 
    122123          pres = f_gradPhi2(ind) 
    123           CALL compute_caldyn_solver(tau,mass,theta,pk,geopot,W, m_il,pres, dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W 
     124          CALL compute_caldyn_solver(tau,phis, mass,theta,pk,geopot,W, m_il,pres, dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W 
    124125       END IF 
    125126       u=f_u(ind) 

  • codes/icosagcm/trunk/src/dynamics/caldyn_kernels_base.F90

    r561 r580  
    4949    Rd = kappa*cpp 
    5050 
    51 #ifdef CPP_DYSL 
    52 !#if 0 
     51    IF(dysl_geopot) THEN 
    5352#include "../kernels/compute_geopot.k90" 
    54 #else 
    55  
     53    ELSE 
    5654    ! Pressure is computed first top-down (temporarily stored in pk) 
    5755    ! Then Exner pressure and geopotential are computed bottom-up 
    … …  
    164162    END IF 
    165163 
    166 #endif 
     164    END IF ! dysl 
    167165 
    168166    !ym flush geopot 
    … …  
    314312    CALL trace_start("compute_caldyn_vert_nh") 
    315313 
    316 #ifdef CPP_DYSL 
     314    IF(dysl) THEN 
    317315!$OMP BARRIER 
    318316#include "../kernels/caldyn_vert_NH.k90" 
    319317!$OMP BARRIER 
    320 #else 
     318    ELSE 
    321319#define ETA_DOT(ij) eta_dot(ij,1) 
    322320#define WCOV(ij) wcov(ij,1) 
    … …  
    365363#undef ETA_DOT 
    366364#undef WCOV 
    367 #endif 
    368  
     365 
     366    END IF ! dysl 
    369367    CALL trace_end("compute_caldyn_vert_nh") 
    370368 

  • codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90

    r561 r580  
    99  PRIVATE 
    1010 
    11   REAL(rstd), PARAMETER :: pbot=1e5, Phi_bot=0., rho_bot=1e6 ! FIXME 
     11  REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6 
    1212 
    1313  LOGICAL, SAVE :: debug_hevi_solver = .FALSE. 
    14   LOGICAL, PARAMETER :: rigid=.TRUE. 
    1514 
    1615  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Coriolis, & 
    … …  
    6261    CALL trace_start("compute_pvort_only")   
    6362!!! Compute shallow-water potential vorticity 
    64 #ifdef CPP_DYSL 
     63    IF(dysl_pvort_only) THEN 
    6564#include "../kernels/pvort_only.k90" 
    66 #else 
     65    ELSE 
     66 
    6767    radius_m2=radius**(-2) 
    6868    DO l = ll_begin,ll_end 
    … …  
    9494 
    9595    ENDDO 
    96 #endif 
     96    
     97    END IF ! dysl 
    9798    CALL trace_end("compute_pvort_only") 
    9899 
    99100  END SUBROUTINE compute_pvort_only 
    100101 
    101   SUBROUTINE compute_NH_geopot(tau, m_ik, m_il, theta, W_il, Phi_il) 
     102  SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il) 
    102103    REAL(rstd),INTENT(IN)    :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs 
     104    REAL(rstd),INTENT(IN)    :: phis(iim*jjm) 
    103105    REAL(rstd),INTENT(IN)    :: m_ik(iim*jjm,llm) 
    104106    REAL(rstd),INTENT(IN)    :: m_il(iim*jjm,llm+1) 
    … …  
    124126    ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1 
    125127 
    126 #ifdef CPP_DYSL 
     128    IF(dysl) THEN 
     129#define PHI_BOT(ij) phis(ij) 
    127130#include "../kernels/compute_NH_geopot.k90" 
    128 #else 
     131#undef PHI_BOT 
     132    ELSE 
    129133!    FIXME : vertical OpenMP parallelism will not work 
    130134    
    … …  
    164168          B_il(ij,1) = A_ik(ij,1) + ml_g2 + tau2_g*rho_bot 
    165169          R_il(ij,1) = ml_g2*( Phi_il(ij,1)-Phi_star_il(ij,1))  & 
    166                + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-Phi_bot) ) 
     170               + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-phis(ij)) ) 
    167171       ENDDO 
    168172       ! inner interfaces 
    … …  
    260264    END DO ! Newton-Raphson 
    261265 
    262 #endif 
     266    END IF ! dysl 
    263267     
    264268  END SUBROUTINE compute_NH_geopot 
    265269 
    266   SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du) 
     270  SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du) 
    267271    REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs 
     272    REAL(rstd),INTENT(IN)    :: phis(iim*jjm) 
    268273    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm) 
    269274    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) 
    … …  
    278283 
    279284    REAL(rstd) :: berni(iim*jjm,llm)         ! (W/m_il)^2 
     285    REAL(rstd) :: berni1(iim*jjm)         ! (W/m_il)^2 
    280286    REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd 
    281287    INTEGER    :: ij, l 
    … …  
    285291    Rd=cpp*kappa 
    286292 
    287 #ifdef CPP_DYSL 
     293    IF(dysl) THEN 
     294 
    288295!$OMP BARRIER 
     296#define PHI_BOT(ij) phis(ij) 
     297#define PHI_BOT_VAR phis 
    289298#include "../kernels/caldyn_solver.k90" 
     299#undef PHI_BOT_VAR 
     300#undef PHI_BOT 
    290301!$OMP BARRIER 
    291 #else 
    292 #define BERNI(ij) berni(ij,1) 
     302 
     303    ELSE 
     304 
     305#define BERNI(ij) berni1(ij) 
    293306    ! FIXME : vertical OpenMP parallelism will not work 
    294307 
    … …  
    310323 
    311324    IF(tau>0) THEN ! solve implicit problem for geopotential 
    312        CALL compute_NH_geopot(tau, rhodz, m_il, theta, W, geopot) 
     325       CALL compute_NH_geopot(tau, phis, rhodz, m_il, theta, W, geopot) 
    313326    END IF 
    314327 
    … …  
    372385    ENDDO 
    373386#undef BERNI 
    374 #endif 
     387 
     388    END IF ! dysl 
    375389 
    376390    CALL trace_end("compute_caldyn_solver") 
    … …  
    396410    Rd=cpp*kappa 
    397411 
    398 #ifdef CPP_DYSL 
     412    IF(dysl_caldyn_fast) THEN 
    399413#include "../kernels/caldyn_fast.k90" 
    400 #else 
     414    ELSE 
     415 
    401416    ! Compute Bernoulli term 
    402417    IF(boussinesq) THEN 
    … …  
    496511       END IF 
    497512    END DO 
    498 #endif        
     513 
     514    END IF ! dysl 
    499515    CALL trace_end("compute_caldyn_fast") 
    500516 
    … …  
    515531    CALL trace_start("compute_caldyn_Coriolis") 
    516532 
    517 #ifdef CPP_DYSL 
     533    IF(dysl_caldyn_coriolis) THEN 
     534 
    518535#include "../kernels/coriolis.k90" 
    519 #else 
     536 
     537    ELSE 
    520538#define FTHETA(ij) Ftheta(ij,1) 
    521539 
    … …  
    653671    END SELECT 
    654672#undef FTHETA 
    655 #endif 
     673 
     674    END IF ! dysl 
    656675 
    657676    CALL trace_end("compute_caldyn_Coriolis") 
    … …  
    667686     
    668687    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
     688    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function 
    669689    REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu 
    670690    INTEGER :: ij,l 
    … …  
    672692    CALL trace_start("compute_caldyn_slow_hydro") 
    673693 
    674 #ifdef CPP_DYSL 
     694    IF(dysl_slow_hydro) THEN 
     695 
    675696#define BERNI(ij,l) berni(ij,l) 
    676697#include "../kernels/caldyn_slow_hydro.k90" 
    677698#undef BERNI 
    678 #else 
    679 #define BERNI(ij) berni(ij,1) 
     699 
     700     ELSE 
     701 
     702#define BERNI(ij) berni1(ij) 
    680703 
    681704    DO l = ll_begin, ll_end 
    … …  
    725748       END IF 
    726749    END DO 
     750 
    727751#undef BERNI 
    728 #endif 
     752 
     753    END IF ! dysl 
    729754    CALL trace_end("compute_caldyn_slow_hydro")     
    730755  END SUBROUTINE compute_caldyn_slow_hydro 
    … …  
    749774    REAL(rstd) :: W_el, W2_el, uu_right, uu_lup, uu_ldown, gPhi2, dP, divG, u2, uu 
    750775 
    751 #ifdef CPP_DYSL 
    752776    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
    753777    REAL(rstd) :: G_el(3*iim*jjm,llm+1) ! horizontal flux of W 
    754778    REAL(rstd) :: v_el(3*iim*jjm,llm+1) 
    755 #else 
    756     REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function 
    757     REAL(rstd) :: G_el(3*iim*jjm) ! horizontal flux of W 
    758     REAL(rstd) :: v_el(3*iim*jjm) 
    759 #endif 
     779 
     780    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function 
     781    REAL(rstd) :: G_el1(3*iim*jjm) ! horizontal flux of W 
     782    REAL(rstd) :: v_el1(3*iim*jjm) 
    760783 
    761784    CALL trace_start("compute_caldyn_slow_NH") 
    762785 
    763 #ifdef CPP_DYSL 
     786    IF(dysl) THEN 
     787 
    764788!$OMP BARRIER 
    765789#include "../kernels/caldyn_slow_NH.k90" 
    766790!$OMP BARRIER 
    767 #else 
     791      
     792     ELSE 
     793 
     794#define BERNI(ij) berni1(ij) 
     795#define G_EL(ij) G_el1(ij) 
     796#define V_EL(ij) v_el1(ij) 
     797 
    768798    DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces) 
    769799       IF(l==1) THEN 
    … …  
    792822          W2_el = .5*le_de(ij+u_right) * & 
    793823               ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) ) 
    794           v_el(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked 
    795           G_el(ij+u_right) = v_el(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el 
     824          V_EL(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked 
     825          G_EL(ij+u_right) = V_EL(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el 
    796826          ! Compute on edge 'lup' 
    797827          W_el  = .5*( W(ij,l)+W(ij+t_lup,l) ) 
    … …  
    800830          W2_el = .5*le_de(ij+u_lup) * & 
    801831               ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) ) 
    802           v_el(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked 
    803           G_el(ij+u_lup) = v_el(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el 
     832          V_EL(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked 
     833          G_EL(ij+u_lup) = V_EL(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el 
    804834          ! Compute on edge 'ldown' 
    805835          W_el  = .5*( W(ij,l)+W(ij+t_ldown,l) ) 
    … …  
    808838          W2_el = .5*le_de(ij+u_ldown) * & 
    809839               ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) ) 
    810           v_el(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked 
    811           G_el(ij+u_ldown) = v_el(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el 
     840          V_EL(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked 
     841          G_EL(ij+u_ldown) = V_EL(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el 
    812842       END DO 
    813843       ! compute GradPhi2, dPhi, dW 
    … …  
    820850               le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 +  & 
    821851               le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 ) 
    822 !          gradPhi2(ij,l) = 0. ! FIXME !! 
    823852 
    824853          dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) -1/(2*Ai(ij))* &  
    825                ( DePhil(ij+u_right,l)*v_el(ij+u_right) + & ! -v.gradPhi,  
    826                  DePhil(ij+u_rup,l)*v_el(ij+u_rup) +     & ! v_el already has le_de 
    827                  DePhil(ij+u_lup,l)*v_el(ij+u_lup) +     & 
    828                  DePhil(ij+u_left,l)*v_el(ij+u_left) +   & 
    829                  DePhil(ij+u_ldown,l)*v_el(ij+u_ldown) + & 
    830                  DePhil(ij+u_rdown,l)*v_el(ij+u_rdown) ) 
     854               ( DePhil(ij+u_right,l)*V_EL(ij+u_right) + & ! -v.gradPhi,  
     855                 DePhil(ij+u_rup,l)*V_EL(ij+u_rup) +     & ! v_el already has le_de 
     856                 DePhil(ij+u_lup,l)*V_EL(ij+u_lup) +     & 
     857                 DePhil(ij+u_left,l)*V_EL(ij+u_left) +   & 
     858                 DePhil(ij+u_ldown,l)*V_EL(ij+u_ldown) + & 
     859                 DePhil(ij+u_rdown,l)*V_EL(ij+u_rdown) ) 
    831860 
    832861          dW(ij,l) = -1./Ai(ij)*(           & ! -div(G_el),  
    833                ne_right*G_el(ij+u_right) +  & ! G_el already has le_de 
    834                ne_rup*G_el(ij+u_rup) +      & 
    835                ne_lup*G_el(ij+u_lup) +      &   
    836                ne_left*G_el(ij+u_left) +    & 
    837                ne_ldown*G_el(ij+u_ldown) +  & 
    838                ne_rdown*G_el(ij+u_rdown)) 
     862               ne_right*G_EL(ij+u_right) +  & ! G_el already has le_de 
     863               ne_rup*G_EL(ij+u_rup) +      & 
     864               ne_lup*G_EL(ij+u_lup) +      &   
     865               ne_left*G_EL(ij+u_left) +    & 
     866               ne_ldown*G_EL(ij+u_ldown) +  & 
     867               ne_rdown*G_EL(ij+u_rdown)) 
    839868       END DO 
    840869    END DO 
    … …  
    843872       ! Compute berni at scalar points 
    844873       DO ij=ij_begin_ext, ij_end_ext 
    845           berni(ij) = & 
     874          BERNI(ij) = & 
    846875               1/(4*Ai(ij))*( & 
    847876                    le_de(ij+u_right)*u(ij+u_right,l)**2 +    & 
    … …  
    860889                    -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1)) 
    861890          hflux(ij+u_right,l) = uu_right*le_de(ij+u_right) 
    862           du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right)) 
     891          du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right)) 
    863892          ! Compute on edge 'lup' 
    864893          uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) & 
    865894                  -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1)) 
    866895          hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup) 
    867           du(ij+u_lup,l) = ne_lup*(berni(ij)-berni(ij+t_lup)) 
     896          du(ij+u_lup,l) = ne_lup*(BERNI(ij)-BERNI(ij+t_lup)) 
    868897          ! Compute on edge 'ldown' 
    869898          uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) & 
    870899                    -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1)) 
    871900          hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown) 
    872           du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown)) 
     901          du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown)) 
    873902       END DO 
    874903    END DO 
    875 #endif 
     904 
     905#undef V_EL 
     906#undef G_EL 
     907#undef BERNI 
     908 
     909    END IF ! dysl 
    876910 
    877911    CALL trace_end("compute_caldyn_slow_NH") 
Note: See TracChangeset for help on using the changeset viewer.