Ignore:
Timestamp:
09/18/17 17:03:05 (7 years ago)
Author:
dubos
Message:

trunk : backported commits 535-540 from devel

Location:
codes/icosagcm/trunk/src/dynamics
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r548 r555  
    3535 
    3636    INTEGER :: i,j,ij,l 
    37     REAL(rstd) :: Rd, p_ik, exner_ik, temp_ik, qv, chi, Rmix 
     37    REAL(rstd) :: Rd, p_ik, exner_ik, temp_ik, qv, chi, Rmix, gv 
    3838    INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext 
    3939 
    … …  
    4747 
    4848    Rd = kappa*cpp 
     49 
     50#ifdef CPP_DYSL 
     51!#if 0 
     52#include "../kernels/compute_geopot.k90" 
     53#else 
    4954 
    5055    ! Pressure is computed first top-down (temporarily stored in pk) 
    … …  
    158163    END IF 
    159164 
     165#endif 
     166 
    160167    !ym flush geopot 
    161168    !$OMP BARRIER 
    … …  
    297304    REAL(rstd),INTENT(INOUT) :: dW(iim*jjm,llm+1) 
    298305    ! local arrays 
    299     REAL(rstd) :: eta_dot(iim*jjm) ! eta_dot in full layers 
    300     REAL(rstd) :: wcov(iim*jjm) ! covariant vertical momentum in full layers 
     306    REAL(rstd) :: eta_dot(iim*jjm, llm) ! eta_dot in full layers 
     307    REAL(rstd) :: wcov(iim*jjm,llm) ! covariant vertical momentum in full layers 
    301308    REAL(rstd) :: W_etadot(iim*jjm,llm) ! vertical flux of vertical momentum 
    302309    ! indices and temporary values 
    … …  
    306313    CALL trace_start("compute_caldyn_vert_nh") 
    307314 
     315#ifdef CPP_DYSL 
     316!#if 0 
     317#include "../kernels/caldyn_vert_NH.k90" 
     318#else 
     319#define ETA_DOT(ij) eta_dot(ij,1) 
     320#define WCOV(ij) wcov(ij,1) 
     321 
    308322    DO l=ll_begin,ll_end 
    309323       ! compute the local arrays 
    … …  
    313327          w_ij = .5*(W(ij,l)+W(ij,l+1))/mass(ij,l) 
    314328          W_etadot(ij,l) = wflux_ij*w_ij 
    315           eta_dot(ij) = wflux_ij / mass(ij,l) 
    316           wcov(ij) = w_ij*(geopot(ij,l+1)-geopot(ij,l)) 
     329          ETA_DOT(ij) = wflux_ij / mass(ij,l) 
     330          WCOV(ij) = w_ij*(geopot(ij,l+1)-geopot(ij,l)) 
    317331       ENDDO 
    318332       ! add NH term to du 
    … …  
    320334      DO ij=ij_begin,ij_end 
    321335          du(ij+u_right,l) = du(ij+u_right,l) & 
    322                - .5*(wcov(ij+t_right)+wcov(ij)) & 
    323                *ne_right*(eta_dot(ij+t_right)-eta_dot(ij)) 
     336               - .5*(WCOV(ij+t_right)+WCOV(ij)) & 
     337               *ne_right*(ETA_DOT(ij+t_right)-ETA_DOT(ij)) 
    324338          du(ij+u_lup,l) = du(ij+u_lup,l) & 
    325                - .5*(wcov(ij+t_lup)+wcov(ij)) & 
    326                *ne_lup*(eta_dot(ij+t_lup)-eta_dot(ij)) 
     339               - .5*(WCOV(ij+t_lup)+WCOV(ij)) & 
     340               *ne_lup*(ETA_DOT(ij+t_lup)-ETA_DOT(ij)) 
    327341          du(ij+u_ldown,l) = du(ij+u_ldown,l) & 
    328                - .5*(wcov(ij+t_ldown)+wcov(ij)) & 
    329                *ne_ldown*(eta_dot(ij+t_ldown)-eta_dot(ij)) 
     342               - .5*(WCOV(ij+t_ldown)+WCOV(ij)) & 
     343               *ne_ldown*(ETA_DOT(ij+t_ldown)-ETA_DOT(ij)) 
    330344       END DO 
    331345    ENDDO 
    … …  
    346360       END DO 
    347361    END DO 
     362 
     363#undef ETA_DOT 
     364#undef WCOV 
     365#endif 
     366 
    348367    CALL trace_end("compute_caldyn_vert_nh") 
    349368 

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

    r548 r555  
    1111  REAL(rstd), PARAMETER :: pbot=1e5, Phi_bot=0., rho_bot=1e6 ! FIXME 
    1212 
    13   LOGICAL, PARAMETER :: debug_hevi_solver = .FALSE. 
     13  LOGICAL, SAVE :: debug_hevi_solver = .FALSE. 
     14  LOGICAL, PARAMETER :: rigid=.TRUE. 
    1415 
    1516  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Coriolis, & 
    … …  
    6162    CALL trace_start("compute_pvort_only")   
    6263!!! Compute shallow-water potential vorticity 
     64#ifdef CPP_DYSL 
     65#include "../kernels/pvort_only.k90" 
     66#else 
    6367    radius_m2=radius**(-2) 
    6468    DO l = ll_begin,ll_end 
    … …  
    9094 
    9195    ENDDO 
    92  
     96#endif 
    9397    CALL trace_end("compute_pvort_only") 
    9498 
    … …  
    114118    REAL(rstd) :: wil, tau2_g, g2, gm2, ml_g2, c2_mik 
    115119 
    116     INTEGER    :: iter, ij, l 
     120    INTEGER    :: iter, ij, l, ij_omp_begin_ext, ij_omp_end_ext 
     121 
     122    CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext) 
     123    ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1 
     124    ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1 
     125 
     126#ifdef CPP_DYSL 
     127!#if 0 
     128#include "../kernels/compute_NH_geopot.k90" 
     129#else 
    117130 
    118131!    FIXME : vertical OpenMP parallelism will not work 
    … …  
    248261 
    249262    END DO ! Newton-Raphson 
     263 
     264#endif 
    250265     
    251266  END SUBROUTINE compute_NH_geopot 
    … …  
    254269    REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs 
    255270    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm) 
    256     REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm) 
     271    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) 
    257272    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm) 
    258273    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) 
    … …  
    263278 
    264279    REAL(rstd) :: m_il(iim*jjm,llm+1)        ! rhodz averaged to interfaces 
    265     REAL(rstd) :: berni(iim*jjm)             ! (W/m_il)^2 
    266     REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij 
     280    REAL(rstd) :: pres(iim*jjm,llm)          ! pressure 
     281    REAL(rstd) :: berni(iim*jjm,llm)             ! (W/m_il)^2 
     282    REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd 
    267283    INTEGER    :: ij, l 
    268284 
    269285    CALL trace_start("compute_caldyn_solver") 
    270286 
     287    Rd=cpp*kappa 
     288 
     289#ifdef CPP_DYSL 
     290!#if 0 
     291#include "../kernels/caldyn_solver.k90" 
     292#else 
     293#define BERNI(ij) berni(ij,1) 
    271294    ! FIXME : vertical OpenMP parallelism will not work 
    272295 
    … …  
    300323       DO ij=ij_begin_ext,ij_end_ext 
    301324          rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l)) 
    302           X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij 
     325          X_ij = (cpp/preff)*kappa*theta(ij,l,1)*rho_ij 
    303326          ! kappa.theta.rho = p/exner  
    304327          ! => X = (p/p0)/(exner/Cp) 
    … …  
    339362       DO ij=ij_begin_ext,ij_end_ext 
    340363          pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow 
    341           berni(ij) = (-.25*g*g)*(              & 
     364          BERNI(ij) = (-.25*g*g)*(              & 
    342365                 (W(ij,l)/m_il(ij,l))**2       & 
    343366               + (W(ij,l+1)/m_il(ij,l+1))**2 ) 
    344367       ENDDO 
    345368       DO ij=ij_begin,ij_end          
    346           du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right)) 
    347           du(ij+u_lup,l)   = ne_lup  *(berni(ij)-berni(ij+t_lup)) 
    348           du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown)) 
     369          du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right)) 
     370          du(ij+u_lup,l)   = ne_lup  *(BERNI(ij)-BERNI(ij+t_lup)) 
     371          du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown)) 
    349372       ENDDO 
    350373    ENDDO 
    351      
     374#undef BERNI 
     375#endif 
     376 
    352377    CALL trace_end("compute_caldyn_solver") 
    353378     
    … …  
    366391 
    367392    INTEGER :: i,j,ij,l 
    368     REAL(rstd) :: Rd, qv, temp, chi, nu, due_right, due_lup, due_ldown 
     393    REAL(rstd) :: Rd, qv, temp, chi, nu, due, due_right, due_lup, due_ldown 
    369394 
    370395    CALL trace_start("compute_caldyn_fast") 
    … …  
    372397    Rd=cpp*kappa 
    373398 
     399#ifdef CPP_DYSL 
     400#include "../kernels/caldyn_fast.k90" 
     401#else 
    374402    ! Compute Bernoulli term 
    375403    IF(boussinesq) THEN 
    … …  
    469497       END IF 
    470498    END DO 
    471         
     499#endif        
    472500    CALL trace_end("compute_caldyn_fast") 
    473501 
    … …  
    482510    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm) 
    483511     
    484     REAL(rstd) :: Ftheta(3*iim*jjm)  ! potential temperature flux 
    485     REAL(rstd) :: uu_right, uu_lup, uu_ldown 
     512    REAL(rstd) :: Ftheta(3*iim*jjm,llm)  ! potential temperature flux 
     513    REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF 
    486514    INTEGER :: ij,iq,l,kdown 
    487515 
    488516    CALL trace_start("compute_caldyn_Coriolis") 
     517 
     518#ifdef CPP_DYSL 
     519!#if 0 
     520#include "../kernels/coriolis.k90" 
     521#else 
     522#define FTHETA(ij) Ftheta(ij,1) 
    489523 
    490524    DO l=ll_begin, ll_end 
    … …  
    493527       !DIR$ SIMD 
    494528          DO ij=ij_begin_ext,ij_end_ext       
    495              Ftheta(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) & 
     529             FTHETA(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) & 
    496530                                  * hflux(ij+u_right,l) 
    497              Ftheta(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) & 
     531             FTHETA(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) & 
    498532                  * hflux(ij+u_lup,l) 
    499              Ftheta(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) & 
     533             FTHETA(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) & 
    500534                  * hflux(ij+u_ldown,l) 
    501535          END DO 
    … …  
    505539             ! dtheta_rhodz =  -div(flux.theta) 
    506540             dtheta_rhodz(ij,l,iq)= & 
    507                   -1./Ai(ij)*(ne_right*Ftheta(ij+u_right)    +  & 
    508                   ne_rup*Ftheta(ij+u_rup)        +  &   
    509                   ne_lup*Ftheta(ij+u_lup)        +  &   
    510                   ne_left*Ftheta(ij+u_left)      +  &   
    511                   ne_ldown*Ftheta(ij+u_ldown)    +  & 
    512                   ne_rdown*Ftheta(ij+u_rdown) ) 
     541                  -1./Ai(ij)*(ne_right*FTHETA(ij+u_right)    +  & 
     542                  ne_rup*FTHETA(ij+u_rup)        +  &   
     543                  ne_lup*FTHETA(ij+u_lup)        +  &   
     544                  ne_left*FTHETA(ij+u_left)      +  &   
     545                  ne_ldown*FTHETA(ij+u_ldown)    +  & 
     546                  ne_rdown*FTHETA(ij+u_rdown) ) 
    513547          END DO 
    514548       END DO 
    … …  
    620654       STOP 
    621655    END SELECT 
     656#undef FTHETA 
     657#endif 
    622658 
    623659    CALL trace_end("compute_caldyn_Coriolis") 
    … …  
    632668    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm) 
    633669     
    634     REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function 
    635     REAL(rstd) :: uu_right, uu_lup, uu_ldown 
     670    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
     671    REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu 
    636672    INTEGER :: ij,l 
    637673 
    638674    CALL trace_start("compute_caldyn_slow_hydro") 
    639675 
    640     le_de(:) = le(:)/de(:) ! FIXME - make sure le_de is what we expect 
     676#ifdef CPP_DYSL 
     677!#if 0 
     678#define BERNI(ij,l) berni(ij,l) 
     679#include "../kernels/caldyn_slow_hydro.k90" 
     680#undef BERNI 
     681#else 
     682#define BERNI(ij) berni(ij,1) 
    641683 
    642684    DO l = ll_begin, ll_end 
    … …  
    658700       !DIR$ SIMD 
    659701       DO ij=ij_begin,ij_end          
    660           berni(ij) = & 
     702          BERNI(ij) = & 
    661703               1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 +    & 
    662704               le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
    … …  
    670712          !DIR$ SIMD 
    671713          DO ij=ij_begin,ij_end 
    672              du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right)) 
    673              du(ij+u_lup,l)   = ne_lup*(berni(ij)-berni(ij+t_lup)) 
    674              du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown)) 
     714             du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right)) 
     715             du(ij+u_lup,l)   = ne_lup*(BERNI(ij)-BERNI(ij+t_lup)) 
     716             du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown)) 
    675717          END DO 
    676718       ELSE 
    … …  
    678720          DO ij=ij_begin,ij_end 
    679721             du(ij+u_right,l) = du(ij+u_right,l) + & 
    680                   ne_right*(berni(ij)-berni(ij+t_right)) 
     722                  ne_right*(BERNI(ij)-BERNI(ij+t_right)) 
    681723             du(ij+u_lup,l)   = du(ij+u_lup,l) + & 
    682                   ne_lup*(berni(ij)-berni(ij+t_lup)) 
     724                  ne_lup*(BERNI(ij)-BERNI(ij+t_lup)) 
    683725             du(ij+u_ldown,l) = du(ij+u_ldown,l) + & 
    684                   ne_ldown*(berni(ij)-berni(ij+t_ldown)) 
     726                  ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown)) 
    685727          END DO 
    686728       END IF 
    687729    END DO 
    688  
     730#undef BERNI 
     731#endif 
    689732    CALL trace_end("compute_caldyn_slow_hydro")     
    690733  END SUBROUTINE compute_caldyn_slow_hydro 
    … …  
    705748    REAL(rstd) :: GradPhi2(3*iim*jjm,llm+1) ! grad_Phi**2 
    706749    REAL(rstd) :: DePhil(3*iim*jjm,llm+1) ! grad(Phi) 
     750     
     751    INTEGER :: ij,l,kdown,kup 
     752    REAL(rstd) :: W_el, W2_el, uu_right, uu_lup, uu_ldown, gPhi2, dP, divG, u2, uu 
     753 
     754#ifdef CPP_DYSL 
     755    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
     756    REAL(rstd) :: G_el(3*iim*jjm,llm+1) ! horizontal flux of W 
     757    REAL(rstd) :: v_el(3*iim*jjm,llm+1) 
     758#else 
    707759    REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function 
    708760    REAL(rstd) :: G_el(3*iim*jjm) ! horizontal flux of W 
    709761    REAL(rstd) :: v_el(3*iim*jjm) 
    710      
    711     INTEGER :: ij,l,kdown,kup 
    712     REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, W2_el 
     762#endif 
    713763 
    714764    CALL trace_start("compute_caldyn_slow_NH") 
    715765 
    716     le_de(:) = le(:)/de(:) ! FIXME - make sure le_de is what we expect 
    717  
     766#ifdef CPP_DYSL 
     767#include "../kernels/caldyn_slow_NH.k90" 
     768#else 
    718769    DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces) 
    719770       IF(l==1) THEN 
    … …  
    789840       END DO 
    790841    END DO 
    791     ! FIXME !! 
    792  !   F_el(:,:)=0.  
    793  !   dPhi(:,:)=0. 
    794  !   dW(:,:)=0. 
    795842 
    796843    DO l=ll_begin, ll_end ! compute on k levels (layers) 
    … …  
    827874       END DO 
    828875    END DO 
    829     ! FIXME !! 
    830     ! du(:,:)=0. 
    831     ! hflux(:,:)=0. 
     876#endif 
    832877 
    833878    CALL trace_end("compute_caldyn_slow_NH") 
Note: See TracChangeset for help on using the changeset viewer.