Changeset 369 for codes


Ignore:
Timestamp:
01/22/16 17:09:54 (8 years ago)
Author:
dubos
Message:

More NH terms - tested with DCMIP3

Location:
codes/icosagcm/trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r366 r369  
    110110       pk = f_pk(ind) 
    111111       geopot = f_geopot(ind) 
     112       du=f_du_fast(ind) 
    112113       IF(hydrostatic) THEN 
     114          du(:,:)=0. 
    113115          CALL compute_geopot(ps,mass,theta, pk,geopot) 
    114116       ELSE 
    … …  
    116118          dW = f_dW_fast(ind) 
    117119          dPhi = f_dPhi_fast(ind) 
    118           CALL compute_caldyn_solver(tau,mass,theta,pk,geopot,W,dPhi,dW) ! computes d(Phi,W)_fast and updates Phi,W 
     120          CALL compute_caldyn_solver(tau,mass,theta,pk,geopot,W,dPhi,dW,du) ! computes d(Phi,W)_fast and updates Phi,W 
    119121       END IF 
    120122       u=f_u(ind) 
    121        du=f_du_fast(ind) 
    122123       CALL compute_caldyn_fast(tau,u,mass,theta,pk,geopot,du) ! computes du_fast and updates du 
    123124    ENDDO 
    … …  
    138139     
    139140    CALL send_message(f_qu,req_qu) ! COM03 
    140      
     141    CALL wait_message(req_qu) ! COM03 
     142        
    141143    DO ind=1,ndomain 
    142144       IF (.NOT. assigned_domain(ind)) CYCLE 
    … …  
    151153       dtheta_rhodz=f_dtheta_rhodz(ind) 
    152154       du=f_du_slow(ind) 
    153        CALL compute_caldyn_slow(u,mass,qu,theta, hflux,convm,dtheta_rhodz,du) ! COM03 
     155       IF(hydrostatic) THEN 
     156          CALL compute_caldyn_slow_hydro(u,mass,hflux,du) 
     157       ELSE 
     158          W = f_W(ind) 
     159          dW = f_dW_slow(ind) 
     160          geopot = f_geopot(ind) 
     161          dPhi = f_dPhi_slow(ind) 
     162          CALL compute_caldyn_slow_NH(u,mass,geopot,W, hflux,du,dPhi,dW) 
     163       END IF 
     164       CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du) 
    154165       IF(caldyn_eta==eta_mass) THEN 
    155166          wflux=f_wflux(ind) 
    156167          wwuu=f_wwuu(ind) 
    157168          dps=f_dps(ind) 
    158           CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du) 
     169          IF(hydrostatic) THEN 
     170             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du) 
     171          ELSE 
     172             PRINT *, 'NH dynamics limited to Lagrangian vertical coordinate so far !'  
     173             STOP 
     174          END IF 
    159175       END IF 
    160176    ENDDO 

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

    r368 r369  
    11MODULE caldyn_kernels_hevi_mod 
    22  USE icosa 
     3  USE trace 
     4  USE omp_para 
     5  USE disvert_mod 
    36  USE transfert_mod 
    47  USE caldyn_kernels_base_mod 
    … …  
    1013  LOGICAL, PARAMETER :: debug_hevi_solver = .FALSE. 
    1114 
    12   PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_slow, & 
     15  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Coriolis, & 
     16       compute_caldyn_slow_hydro, compute_caldyn_slow_NH, & 
    1317       compute_caldyn_solver, compute_caldyn_fast 
    1418 
    … …  
    1620 
    1721  SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta) 
    18     USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop 
    19     USE exner_mod 
    20     USE trace 
    21     USE omp_para 
    22     IMPLICIT NONE 
    2322    REAL(rstd),INTENT(IN)  :: ps(iim*jjm) 
    2423    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm) 
    … …  
    5554 
    5655  SUBROUTINE compute_pvort_only(u,rhodz,qu,qv) 
    57     USE exner_mod 
    58     USE trace 
    59     USE omp_para 
    60     IMPLICIT NONE 
    6156    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
    6257    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm) 
    … …  
    124119 
    125120  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 
    129121    REAL(rstd),INTENT(IN)    :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs 
    130122    REAL(rstd),INTENT(IN)    :: m_ik(iim*jjm,llm) 
    … …  
    282274  END SUBROUTINE compute_NH_geopot 
    283275 
    284   SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk, geopot,W, dPhi,dW) 
    285     USE disvert_mod 
    286     USE exner_mod 
    287     USE trace 
    288     USE omp_para 
    289     IMPLICIT NONE 
     276  SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk, geopot,W, dPhi,dW,du) 
    290277    REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs 
    291278    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm) 
    … …  
    296283    REAL(rstd),INTENT(OUT)   :: dW(iim*jjm,llm+1) 
    297284    REAL(rstd),INTENT(OUT)   :: dPhi(iim*jjm,llm+1) 
     285    REAL(rstd),INTENT(OUT)   :: du(3*iim*jjm,llm) 
    298286 
    299287    REAL(rstd) :: m_il(iim*jjm,llm+1)        ! rhodz averaged to interfaces 
     288    REAL(rstd) :: berni(iim*jjm)             ! (W/m_il)^2 
    300289    REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij 
    301290    INTEGER    :: ij, l 
    … …  
    342331    ENDDO 
    343332 
    344     ! Update W, compute tendencies for geopotential and vertical momentum 
     333    ! Update W, compute tendencies 
    345334    DO l=2,llm 
    346335       !DIR$ SIMD 
    … …  
    368357    !    ENDDO 
    369358     
    370     ! Compute Exner function (needed by compute_caldyn_fast) 
     359    ! Compute Exner function (needed by compute_caldyn_fast) and du=g^-2.grad(w^2) 
    371360    DO l=1,llm 
    372361       !DIR$ SIMD 
    373362       DO ij=ij_begin_ext,ij_end_ext 
    374363          pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow 
     364          berni(ij) = (-.5*g*g)*( (W(ij,l)/m_il(ij,l))**2 & 
     365               + (W(ij,l+1)/m_il(ij,l+1))**2 ) 
     366       ENDDO 
     367       DO ij=ij_begin,ij_end          
     368          du(ij+u_right,l) = ne_right*berni(ij)+ne_left *berni(ij+t_right) 
     369          du(ij+u_lup,l)   = ne_lup  *berni(ij)+ne_rdown*berni(ij+t_lup) 
     370          du(ij+u_ldown,l) = ne_ldown*berni(ij)+ne_rup  *berni(ij+t_ldown) 
    375371       ENDDO 
    376372    ENDDO 
    … …  
    381377   
    382378  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du) 
    383     USE icosa 
    384     USE disvert_mod 
    385     USE exner_mod 
    386     USE trace 
    387     USE omp_para 
    388     IMPLICIT NONE 
    389379    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs 
    390380    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! OUT if tau>0 
    … …  
    393383    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) 
    394384    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) 
    395     REAL(rstd),INTENT(OUT)   :: du(iim*3*jjm,llm) 
     385    REAL(rstd),INTENT(INOUT)   :: du(iim*3*jjm,llm) 
    396386    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
    397387 
    … …  
    424414    END IF ! Boussinesq/compressible 
    425415 
    426 !!! u:=u+tau*du, du = gradients of Bernoulli and Exner functions 
     416!!! u:=u+tau*du, du = -grad(B)-theta.grad(pi) 
    427417    DO l=ll_begin,ll_end 
    428418       !DIR$ SIMD 
    … …  
    437427                         *(ne_ldown*pk(ij,l)   +ne_rup*pk(ij+t_ldown,l))      & 
    438428                        +  ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) 
    439           IF(.NOT.DEC) THEN 
    440              due_right = due_right/de(ij+u_right) 
    441              due_lup   = due_lup/de(ij+u_lup) 
    442              due_ldown = due_ldown/de(ij+u_ldown) 
    443           END IF 
    444           du(ij+u_right,l) = due_right 
    445           du(ij+u_lup,l)   = due_lup 
    446           du(ij+u_ldown,l) = due_ldown 
    447           u(ij+u_right,l) = u(ij+u_right,l) + tau*due_right 
    448           u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*due_lup 
    449           u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*due_ldown 
     429          du(ij+u_right,l) = du(ij+u_right,l) + due_right 
     430          du(ij+u_lup,l)   = du(ij+u_lup,l)   + due_lup 
     431          du(ij+u_ldown,l) = du(ij+u_ldown,l) + due_ldown 
     432          u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l) 
     433          u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l) 
     434          u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l) 
    450435       ENDDO 
    451436    ENDDO 
    … …  
    455440  END SUBROUTINE compute_caldyn_fast 
    456441 
    457   SUBROUTINE compute_caldyn_slow(u,rhodz,qu,theta, hflux,convm, dtheta_rhodz, du) 
    458     USE icosa 
    459     USE disvert_mod 
    460     USE exner_mod 
    461     USE trace 
    462     USE omp_para 
    463     IMPLICIT NONE 
    464     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)    ! prognostic "velocity" 
    465     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    466     REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm) 
     442  SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du) 
     443    REAL(rstd),INTENT(IN)  :: hflux(3*iim*jjm,llm)  ! hflux in kg/s 
    467444    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature 
    468  
    469     REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s 
     445    REAL(rstd),INTENT(IN)  :: qu(3*iim*jjm,llm) 
    470446    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence 
    471447    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) 
    472     REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) 
    473  
    474     REAL(rstd) :: cor_NT(iim*jjm,llm)  ! NT coriolis force u.(du/dPhi) 
    475     REAL(rstd) :: urel(3*iim*jjm,llm)  ! relative velocity 
    476     REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux 
    477     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
    478  
    479     INTEGER :: ij,l 
     448    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm) 
     449     
     450    REAL(rstd) :: Ftheta(3*iim*jjm)  ! potential temperature flux 
    480451    REAL(rstd) :: uu_right, uu_lup, uu_ldown 
    481  
    482     CALL trace_start("compute_caldyn_slow") 
    483  
    484     DO l = ll_begin, ll_end 
    485 !!!  Compute mass and theta fluxes 
    486        IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    487        !DIR$ SIMD 
    488        DO ij=ij_begin_ext,ij_end_ext 
    489           uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) 
    490           uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) 
    491           uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) 
    492           IF(DEC) THEN 
    493              uu_right= uu_right*le(ij+u_right)/de(ij+u_right) 
    494              uu_lup  = uu_lup  *le(ij+u_lup)/de(ij+u_lup) 
    495              uu_ldown= uu_ldown*le(ij+u_ldown)/de(ij+u_ldown) 
    496           ELSE 
    497              uu_right= uu_right*le(ij+u_right) 
    498              uu_lup  = uu_lup  *le(ij+u_lup) 
    499              uu_ldown= uu_ldown*le(ij+u_ldown) 
    500           END IF 
    501           hflux(ij+u_right,l)=uu_right 
    502           hflux(ij+u_lup,l)  =uu_lup 
    503           hflux(ij+u_ldown,l)=uu_ldown 
    504           Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*uu_right 
    505           Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*uu_lup 
    506           Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*uu_ldown 
    507        ENDDO 
    508  
    509 !!! compute horizontal divergence of fluxes 
    510        !DIR$ SIMD 
     452    INTEGER :: ij,l,kdown 
     453 
     454    CALL trace_start("compute_caldyn_Coriolis") 
     455 
     456    DO l=ll_begin, ll_end 
     457       ! compute theta flux 
     458       !DIR$ SIMD 
     459       DO ij=ij_begin_ext,ij_end_ext       
     460          Ftheta(ij+u_right) = 0.5*(theta(ij,l)+theta(ij+t_right,l)) & 
     461                                  * hflux(ij+u_right,l) 
     462          Ftheta(ij+u_lup)   = 0.5*(theta(ij,l)+theta(ij+t_lup,l)) & 
     463                                  * hflux(ij+u_lup,l) 
     464          Ftheta(ij+u_ldown) = 0.5*(theta(ij,l)+theta(ij+t_ldown,l)) & 
     465                                  * hflux(ij+u_ldown,l) 
     466       ENDDO 
     467       ! compute horizontal divergence of fluxes 
    511468       DO ij=ij_begin,ij_end          
    512469          ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 
    … …  
    517474               ne_ldown*hflux(ij+u_ldown,l)   +  &   
    518475               ne_rdown*hflux(ij+u_rdown,l))     
    519  
    520           ! signe ? attention d (rho theta dz) 
    521476          ! dtheta_rhodz =  -div(flux.theta) 
    522           dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  & 
    523                ne_rup*Ftheta(ij+u_rup,l)        +  &   
    524                ne_lup*Ftheta(ij+u_lup,l)        +  &   
    525                ne_left*Ftheta(ij+u_left,l)      +  &   
    526                ne_ldown*Ftheta(ij+u_ldown,l)    +  &   
    527                ne_rdown*Ftheta(ij+u_rdown,l)) 
    528        ENDDO 
    529  
     477          dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right)    +  & 
     478               ne_rup*Ftheta(ij+u_rup)        +  &   
     479               ne_lup*Ftheta(ij+u_lup)        +  &   
     480               ne_left*Ftheta(ij+u_left)      +  &   
     481               ne_ldown*Ftheta(ij+u_ldown)    +  &   
     482               ne_rdown*Ftheta(ij+u_rdown)) 
     483       ENDDO 
    530484    END DO 
    531485 
    532486!!! Compute potential vorticity (Coriolis) contribution to du 
    533  
    534487    SELECT CASE(caldyn_conserv) 
     488 
    535489    CASE(energy) ! energy-conserving TRiSK 
    536  
    537        CALL wait_message(req_qu) 
    538490 
    539491       DO l=ll_begin,ll_end 
    … …  
    573525                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      & 
    574526                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l)) 
    575              IF(DEC) THEN 
    576                 du(ij+u_right,l) = .5*uu_right 
    577                 du(ij+u_lup,l)   = .5*uu_lup 
    578                 du(ij+u_ldown,l)   = .5*uu_ldown 
    579              ELSE 
    580                 du(ij+u_right,l) = .5*uu_right/de(ij+u_right) 
    581                 du(ij+u_lup,l)   = .5*uu_lup  /de(ij+u_lup) 
    582                 du(ij+u_ldown,l) = .5*uu_ldown/de(ij+u_ldown) 
    583              END IF 
     527             du(ij+u_right,l) = du(ij+u_right,l) + .5*uu_right 
     528             du(ij+u_lup,l)   = du(ij+u_lup,l)   + .5*uu_lup 
     529             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + .5*uu_ldown 
    584530          ENDDO 
    585531       ENDDO 
    … …  
    623569                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              & 
    624570                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l) 
    625              IF(DEC) THEN 
    626                 du(ij+u_right,l) = qu(ij+u_right,l)*uu_right 
    627                 du(ij+u_lup,l)   = qu(ij+u_lup,l)  *uu_lup 
    628                 du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu_ldown 
    629              ELSE 
    630                 du(ij+u_right,l) = qu(ij+u_right,l)*uu_right/de(ij+u_right) 
    631                 du(ij+u_lup,l)   = qu(ij+u_lup,l)  *uu_lup  /de(ij+u_lup) 
    632                 du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu_ldown/de(ij+u_ldown) 
    633              END IF 
    634           ENDDO 
    635        ENDDO 
    636  
     571 
     572             du(ij+u_right,l) = du(ij+u_right,l) + .5*uu_right 
     573             du(ij+u_lup,l)   = du(ij+u_lup,l)   + .5*uu_lup 
     574             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + .5*uu_ldown 
     575          END DO 
     576       END DO 
    637577    CASE DEFAULT 
    638578       STOP 
    639579    END SELECT 
    640580 
    641     ! Compute bernouilli term = Kinetic Energy 
     581    CALL trace_end("compute_caldyn_Coriolis") 
     582 
     583  END SUBROUTINE compute_caldyn_Coriolis 
     584 
     585  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du) 
     586    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity" 
     587    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
     588    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s 
     589    REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm) 
     590     
     591    REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function 
     592    REAL(rstd) :: uu_right, uu_lup, uu_ldown 
     593    INTEGER :: ij,l 
     594 
     595    CALL trace_start("compute_caldyn_slow_hydro") 
     596 
    642597    le_de(:) = le(:)/de(:) ! FIXME - make sure le_de is what we expect 
    643     DO l=ll_begin,ll_end 
     598 
     599    DO l = ll_begin, ll_end 
     600       !  Compute mass fluxes 
     601       IF (caldyn_conserv==energy) CALL test_message(req_qu)  
     602       !DIR$ SIMD 
     603       DO ij=ij_begin_ext,ij_end_ext 
     604          uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) 
     605          uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) 
     606          uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) 
     607          uu_right= uu_right*le_de(ij+u_right) 
     608          uu_lup  = uu_lup  *le_de(ij+u_lup) 
     609          uu_ldown= uu_ldown*le_de(ij+u_ldown) 
     610          hflux(ij+u_right,l)=uu_right 
     611          hflux(ij+u_lup,l)  =uu_lup 
     612          hflux(ij+u_ldown,l)=uu_ldown 
     613       ENDDO 
     614       ! Compute Bernoulli=kinetic energy 
    644615       !DIR$ SIMD 
    645616       DO ij=ij_begin,ij_end          
    646           IF(DEC) THEN 
    647              berni(ij,l) = & 
    648                   1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 +    & 
    649                   le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
    650                   le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
    651                   le_de(ij+u_left)*u(ij+u_left,l)**2 +       & 
    652                   le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
    653                   le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    654           ELSE 
    655              berni(ij,l) = & 
    656                   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    & 
    657                   le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
    658                   le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
    659                   le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       & 
    660                   le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
    661                   le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    662           END IF 
    663        ENDDO 
    664     ENDDO 
    665  
    666 !!! Add gradients of Bernoulli and Exner functions to du 
    667     DO l=ll_begin,ll_end 
     617          berni(ij) = & 
     618               1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 +    & 
     619               le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
     620               le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
     621               le_de(ij+u_left)*u(ij+u_left,l)**2 +       & 
     622               le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
     623               le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
     624       ENDDO 
     625       ! Compute du=grad(Bernoulli) 
    668626       !DIR$ SIMD 
    669627       DO ij=ij_begin,ij_end 
    670           IF(DEC) THEN 
    671              du(ij+u_right,l) = du(ij+u_right,l)                       & 
    672                   + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) 
    673              du(ij+u_lup,l) = du(ij+u_lup,l)                           & 
    674                   + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) 
    675              du(ij+u_ldown,l) = du(ij+u_ldown,l)                       & 
    676                   + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) 
    677           ELSE 
    678              du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right)            & 
    679                   * ( ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) ) 
    680              du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup)                 & 
    681                   * ( ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) ) 
    682              du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown)            & 
    683                   * ( ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) ) 
    684           END IF 
     628             du(ij+u_right,l) = ne_right*berni(ij)+ne_left*berni(ij+t_right) 
     629             du(ij+u_lup,l)   = ne_lup*berni(ij)+ne_rdown*berni(ij+t_lup) 
     630             du(ij+u_ldown,l) = ne_ldown*berni(ij)+ne_rup*berni(ij+t_ldown) 
    685631       END DO 
    686     ENDDO 
    687  
    688     CALL trace_end("compute_caldyn_slow") 
    689  
    690   END SUBROUTINE compute_caldyn_slow 
     632    END DO 
     633 
     634    CALL trace_end("compute_caldyn_slow_hydro")     
     635  END SUBROUTINE compute_caldyn_slow_hydro 
     636 
     637  SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, hflux,du,dPhi,dW) 
     638    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity" 
     639    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)  ! rho*dz 
     640    REAL(rstd),INTENT(IN)  :: Phi(iim*jjm,llm+1)  ! prognostic geopotential 
     641    REAL(rstd),INTENT(IN)  :: W(iim*jjm,llm+1)    ! prognostic vertical momentum 
     642 
     643    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s 
     644    REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm) 
     645    REAL(rstd),INTENT(OUT) :: dW(iim*jjm,llm+1) 
     646    REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1) 
     647 
     648    REAL(rstd) :: w_il(3*iim*jjm,llm+1) ! Wil/mil 
     649    REAL(rstd) :: F_el(3*iim*jjm,llm+1) ! NH mass flux 
     650    REAL(rstd) :: GradPhi2(3*iim*jjm,llm+1) ! grad_Phi**2 
     651    REAL(rstd) :: DePhil(3*iim*jjm,llm+1) ! grad(Phi) 
     652    REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function 
     653    REAL(rstd) :: G_el(3*iim*jjm) ! horizontal flux of W 
     654    REAL(rstd) :: v_el(3*iim*jjm) 
     655     
     656    INTEGER :: ij,l,kdown,kup 
     657    REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, W2_el 
     658 
     659    CALL trace_start("compute_caldyn_slow_NH") 
     660 
     661    le_de(:) = le(:)/de(:) ! FIXME - make sure le_de is what we expect 
     662 
     663    DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces) 
     664       IF(l==1) THEN 
     665          kdown=1 
     666       ELSE 
     667          kdown=l-1 
     668       END IF 
     669       IF(l==llm+1) THEN 
     670          kup=llm 
     671       ELSE 
     672          kup=l 
     673       END IF 
     674       ! compute mil, wil=Wil/mil 
     675       DO ij=ij_begin_ext, ij_end_ext 
     676          w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup)) 
     677       END DO 
     678       ! compute DePhi, v_el, G_el, F_el 
     679       ! v_el, W2_el and therefore G_el incorporate metric factor le_de 
     680       ! while DePhil, W_el and F_el don't 
     681       DO ij=ij_begin_ext, ij_end_ext 
     682          ! Compute on edge 'right' 
     683          W_el  = .5*( W(ij,l)+W(ij+t_right,l) ) 
     684          DePhil(ij+u_right,l) = ne_right*(Phi(ij+t_right,l)-Phi(ij,l)) 
     685          F_el(ij+u_right,l)   = DePhil(ij+u_right,l)*W_el 
     686          W2_el = .5*le_de(ij+u_right) * & 
     687               ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) ) 
     688          v_el(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) 
     689          G_el(ij+u_right) = v_el(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el 
     690          ! Compute on edge 'lup' 
     691          W_el  = .5*( W(ij,l)+W(ij+t_lup,l) ) 
     692          DePhil(ij+u_lup,l) = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l)) 
     693          F_el(ij+u_lup,l)   = DePhil(ij+u_lup,l)*W_el 
     694          W2_el = .5*le_de(ij+u_lup) * & 
     695               ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) ) 
     696          v_el(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown) ) 
     697          G_el(ij+u_lup) = v_el(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el 
     698          ! Compute on edge 'ldown' 
     699          W_el  = .5*( W(ij,l)+W(ij+t_ldown,l) ) 
     700          DePhil(ij+u_ldown,l) = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l)) 
     701          F_el(ij+u_ldown,l) = DePhil(ij+u_ldown,l)*W_el 
     702          W2_el = .5*le_de(ij+u_ldown) * & 
     703               ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) ) 
     704          v_el(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown) ) 
     705          G_el(ij+u_ldown) = v_el(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el 
     706       END DO 
     707       ! compute GradPhi2, dPhi, dW 
     708       DO ij=ij_begin_ext, ij_end_ext 
     709          gradPhi2(ij,l) = & 
     710               1/(2*Ai(ij))*(le_de(ij+u_right)*DePhil(ij+u_right,l)**2 + & 
     711               le_de(ij+u_rup)*DePhil(ij+u_rup,l)**2 +      & 
     712               le_de(ij+u_lup)*DePhil(ij+u_lup,l)**2 +      & 
     713               le_de(ij+u_left)*DePhil(ij+u_left,l)**2 +    & 
     714               le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 +  & 
     715               le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 ) 
     716          dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l)-1/(2*Ai(ij))* &  
     717               ( DePhil(ij+u_right,l)*v_el(ij+u_right) + & ! -v.gradPhi,  
     718                 DePhil(ij+u_rup,l)*v_el(ij+u_rup) +     & ! v_el already has le_de 
     719                 DePhil(ij+u_lup,l)*v_el(ij+u_lup) +     & 
     720                 DePhil(ij+u_left,l)*v_el(ij+u_left) +   & 
     721                 DePhil(ij+u_ldown,l)*v_el(ij+u_ldown) + & 
     722                 DePhil(ij+u_rdown,l)*v_el(ij+u_rdown) ) 
     723          dW(ij,l) = -1./Ai(ij)*(           & ! -div(G_el),  
     724               ne_right*G_el(ij+u_right) +  & ! G_el already has le_de 
     725               ne_rup*G_el(ij+u_rup) +      & 
     726               ne_lup*G_el(ij+u_lup) +      &   
     727               ne_left*G_el(ij+u_left) +    & 
     728               ne_ldown*G_el(ij+u_ldown) +  & 
     729               ne_rdown*G_el(ij+u_rdown)) 
     730       END DO 
     731    END DO 
     732     
     733    DO l=ll_begin, ll_end ! compute on k levels (layers) 
     734       ! Compute berni at scalar points 
     735       DO ij=ij_begin_ext, ij_end_ext 
     736          berni(ij) = & 
     737               1/(4*Ai(ij))*( & 
     738                    le_de(ij+u_right)*u(ij+u_right,l)**2 +    & 
     739                    le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
     740                    le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
     741                    le_de(ij+u_left)*u(ij+u_left,l)**2 +       & 
     742                    le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
     743                    le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) & 
     744               - .5*( gradPhi2(ij,l)  *w_il(ij,l)**2 +   & 
     745                      gradPhi2(ij,l+1)*w_il(ij,l+1)**2 ) 
     746       END DO 
     747       ! Compute mass flux and grad(berni) at edges 
     748       DO ij=ij_begin_ext, ij_end_ext 
     749          ! Compute on edge 'right' 
     750          uu_right = 0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) & 
     751                    -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1)) 
     752          hflux(ij+u_right,l) = uu_right*le_de(ij+u_right) 
     753          du(ij+u_right,l) = ne_right*berni(ij)+ne_left*berni(ij+t_right) 
     754          ! Compute on edge 'lup' 
     755          uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) & 
     756                  -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1)) 
     757          hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup) 
     758          du(ij+u_lup,l) = ne_lup*berni(ij)+ne_rdown*berni(ij+t_lup) 
     759          ! Compute on edge 'ldown' 
     760          uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) & 
     761                    -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1)) 
     762          hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown) 
     763          du(ij+u_ldown,l) = ne_ldown*berni(ij)+ne_rup*berni(ij+t_ldown) 
     764       END DO 
     765    END DO 
     766 
     767    CALL trace_end("compute_caldyn_slow_NH") 
     768 
     769  END SUBROUTINE compute_caldyn_slow_NH 
    691770 
    692771END MODULE caldyn_kernels_hevi_mod 

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

    r366 r369  
    12301230                                Om      = 0.d0,                 &       ! Rotation Rate of Earth 
    12311231                                as      = a/X,                  &       ! New Radius of small Earth      
    1232 !                               u0      = 20.d0,                &       ! Reference Velocity  
    1233                                 u0      = 0.d0,                 &       ! FIXME : no zonal wind for NH tests 
     1232                                u0      = 20.d0,                &       ! Reference Velocity  
     1233!                               u0      = 0.d0,                 &       ! FIXME : no zonal wind for NH tests 
    12341234                                Teq     = 300.d0,               &       ! Temperature at Equator     
    12351235                                Peq     = 100000.d0,            &       ! Reference PS at Equator 
Note: See TracChangeset for help on using the changeset viewer.