Ignore:
Timestamp:
08/27/18 18:35:24 (6 years ago)
Author:
dubos
Message:

devel : Gassmann (2018) modification of TRiSK Coriolis discretization

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r734 r735  
    5050  END SUBROUTINE compute_theta 
    5151 
    52   SUBROUTINE compute_pvort_only(u,rhodz,qu,qv) 
     52  SUBROUTINE compute_pvort_only(u,rhodz,qu,qv,hv_) 
    5353    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
    5454    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm) 
    5555    REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm) 
    5656    REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm) 
     57    REAL(rstd),INTENT(OUT) :: hv_(iim*2*jjm,llm) 
    5758 
    5859    INTEGER :: ij,l 
    … …  
    7677               + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 
    7778          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv 
     79          hv_(ij+z_up,l) = hv 
    7880           
    7981          etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)   & 
    … …  
    8486               + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 
    8587          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 
     88          hv_(ij+z_down,l) = hv 
    8689       ENDDO 
    8790 
    … …  
    668671       ENDDO 
    669672 
     673    CASE(conserv_gassmann) ! energy-conserving TRiSK modified by Gassmann (2018) 
     674 
     675       DO l=ll_begin,ll_end 
     676          !DIR$ SIMD 
     677          DO ij=ij_begin,ij_end          
     678             uu_right = & 
     679                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)  *qu(ij+t_right+u_lup,l)+                 & 
     680                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)  *qu(ij+u_rup,l)+                         & 
     681                .5*wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+     & 
     682                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*qu(ij+u_rdown,l)+                       & 
     683                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*qu(ij+t_right+u_ldown,l)+               & 
     684                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*qu(ij+u_rdown,l)+               & 
     685                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*qu(ij+t_right+u_ldown,l)+       & 
     686               .5*wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         & 
     687                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*qu(ij+t_right+u_lup,l)+           & 
     688                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*qu(ij+u_rup,l) 
     689             uu_lup = & 
     690                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*qu(ij+t_lup+u_ldown,l) +                   & 
     691                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*qu(ij+u_left,l) +                         & 
     692               .5*wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +       & 
     693                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*qu(ij+u_rup,l) +                          & 
     694                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*qu(ij+t_lup+u_right,l)+                     &  
     695                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*qu(ij+u_rup,l) +                   & 
     696                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*qu(ij+t_lup+u_right,l) +              & 
     697               .5*wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) + & 
     698                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*qu(ij+t_lup+u_ldown,l) +            & 
     699                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*qu(ij+u_left,l) 
     700             uu_ldown = & 
     701                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*qu(ij+t_ldown,l+u_right) +               & 
     702                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*qu(ij+u_rdown,l) +                       & 
     703               .5*wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +        & 
     704                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*qu(ij+u_left,l) +                          & 
     705                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*qu(ij+t_ldown+u_lup,l) +                  &  
     706                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*qu(ij+u_left,l) +                    & 
     707                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*qu(ij+t_ldown+u_lup,l) +          & 
     708               .5*wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      & 
     709                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*qu(ij+t_ldown+u_right,l) +      & 
     710                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*qu(ij+u_rdown,l) 
     711             du(ij+u_right,l) = du(ij+u_right,l) + uu_right 
     712             du(ij+u_lup,l)   = du(ij+u_lup,l)   + uu_lup 
     713             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + uu_ldown 
     714          ENDDO 
     715       ENDDO 
     716 
    670717    CASE(conserv_enstrophy) ! enstrophy-conserving TRiSK 
    671718 
    … …  
    723770  END SUBROUTINE compute_caldyn_Coriolis 
    724771 
    725   SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,Kv,du, zero) 
     772  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hv, hflux,Kv,du, zero) 
    726773    LOGICAL, INTENT(IN) :: zero 
    727774    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity" 
    728775    REAL(rstd),INTENT(IN)  :: Kv(2*iim*jjm,llm)   ! kinetic energy at vertices 
     776    REAL(rstd),INTENT(IN)  :: hv(2*iim*jjm,llm)   ! height/mass averaged to vertices 
    729777    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    730778    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s 
    … …  
    751799       !  Compute mass fluxes 
    752800       IF (caldyn_conserv==conserv_energy) CALL test_message(req_qu)  
    753        !DIR$ SIMD 
    754        DO ij=ij_begin_ext,ij_end_ext 
    755           uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) 
    756           uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) 
    757           uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) 
    758           uu_right= uu_right*le_de(ij+u_right) 
    759           uu_lup  = uu_lup  *le_de(ij+u_lup) 
    760           uu_ldown= uu_ldown*le_de(ij+u_ldown) 
    761           hflux(ij+u_right,l)=uu_right 
    762           hflux(ij+u_lup,l)  =uu_lup 
    763           hflux(ij+u_ldown,l)=uu_ldown 
    764        ENDDO 
     801 
     802       IF(caldyn_kinetic==kinetic_trisk) THEN 
     803          !DIR$ SIMD 
     804          DO ij=ij_begin_ext,ij_end_ext 
     805             uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) 
     806             uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) 
     807             uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) 
     808             uu_right= uu_right*le_de(ij+u_right) 
     809             uu_lup  = uu_lup  *le_de(ij+u_lup) 
     810             uu_ldown= uu_ldown*le_de(ij+u_ldown) 
     811             hflux(ij+u_right,l)=uu_right 
     812             hflux(ij+u_lup,l)  =uu_lup 
     813             hflux(ij+u_ldown,l)=uu_ldown 
     814          ENDDO 
     815       ELSE ! mass flux deriving from consistent kinetic energy 
     816          !DIR$ SIMD 
     817          DO ij=ij_begin_ext,ij_end_ext 
     818             uu_right=0.5*(hv(ij+z_rup,l)+hv(ij+z_rdown,l))*u(ij+u_right,l) 
     819             uu_lup=0.5*(hv(ij+z_up,l)+hv(ij+z_lup,l))*u(ij+u_lup,l) 
     820             uu_ldown=0.5*(hv(ij+z_ldown,l)+hv(ij+z_down,l))*u(ij+u_ldown,l) 
     821             uu_right= uu_right*le_de(ij+u_right) 
     822             uu_lup  = uu_lup  *le_de(ij+u_lup) 
     823             uu_ldown= uu_ldown*le_de(ij+u_ldown) 
     824             hflux(ij+u_right,l)=uu_right 
     825             hflux(ij+u_lup,l)  =uu_lup 
     826             hflux(ij+u_ldown,l)=uu_ldown 
     827          ENDDO 
     828       END IF 
     829 
    765830       ! Compute Bernoulli=kinetic energy 
    766831       IF(caldyn_kinetic==kinetic_trisk) THEN 
Note: See TracChangeset for help on using the changeset viewer.