Changeset 735 for codes


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

devel : Gassmann (2018) modification of TRiSK Coriolis discretization

Location:
codes/icosagcm/devel/src/dynamics
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r733 r735  
    6161    CASE('energy') 
    6262       caldyn_conserv=conserv_energy 
     63    CASE('energy_gassmann') 
     64       caldyn_conserv=conserv_gassmann 
    6365    CASE('enstrophy') 
    6466       caldyn_conserv=conserv_enstrophy 
    6567    CASE DEFAULT 
    6668       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', & 
    67             TRIM(def),'> options are <energy>, <enstrophy>' 
     69            TRIM(def),'> options are <energy>, <energy_gassmann>, <enstrophy>' 
    6870       STOP 
    6971    END SELECT 
    … …  
    156158 
    157159  SUBROUTINE allocate_caldyn 
    158     CALL allocate_field(f_qu,field_u,type_real,llm)  
    159     CALL allocate_field(f_qv,field_z,type_real,llm)  
    160     CALL allocate_field(f_Kv,field_z,type_real,llm)  
     160    CALL allocate_field(f_qu,field_u,type_real,llm, name='qu') 
     161    CALL allocate_field(f_qv,field_z,type_real,llm, name='qv') 
     162    CALL allocate_field(f_Kv,field_z,type_real,llm, name='Kv') 
     163    CALL allocate_field(f_hv,field_z,type_real,llm, name='hv') 
    161164    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk') 
    162     CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')     
     165    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu') 
    163166    CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a 
    164167    IF(.NOT.hydrostatic) THEN 

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

    r734 r735  
    5151    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:) 
    5252    REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:) 
    53     REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:), Kv(:,:) 
     53    REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:), Kv(:,:), hv(:,:) 
    5454 
    5555! temporary shared variable 
    … …  
    140140       qu=f_qu(ind) 
    141141       qv=f_qv(ind) 
    142        CALL compute_pvort_only(u,mass,qu,qv) 
     142       hv=f_hv(ind) 
     143       CALL compute_pvort_only(u,mass,qu,qv,hv) 
    143144       IF(caldyn_kinetic==kinetic_consistent) THEN 
    144145          Kv=f_Kv(ind) 
    … …  
    151152 
    152153    IF(caldyn_kinetic==kinetic_consistent) THEN 
     154       CALL transfert_request(f_hv,req_z1_scal) 
    153155       CALL send_message(f_Kv,req_Kv) 
    154156       CALL wait_message(req_Kv) 
    … …  
    169171 
    170172       IF(hydrostatic) THEN 
     173          hv=f_hv(ind) 
    171174          Kv=f_Kv(ind) 
    172           CALL compute_caldyn_slow_hydro(u,mass,hflux,Kv,du, .TRUE.) 
     175          CALL compute_caldyn_slow_hydro(u,mass,hv, hflux,Kv,du, .TRUE.) 
    173176       ELSE 
    174177          W = f_W(ind) 

  • 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 

  • codes/icosagcm/devel/src/dynamics/caldyn_vars.f90

    r734 r735  
    66  SAVE 
    77 
    8   INTEGER, PARAMETER :: conserv_energy=1, conserv_enstrophy=2, & 
     8  INTEGER, PARAMETER :: conserv_energy=1, conserv_enstrophy=2, conserv_gassmann=3, & 
    99       kinetic_trisk=1, kinetic_consistent=2, & 
    1010       caldyn_vert_noncons=1, caldyn_vert_cons=2, & 
    … …  
    2222 
    2323  ! temporary shared variables for caldyn 
    24   TYPE(t_field),POINTER :: f_qu(:), f_qv(:), f_Kv(:), & 
     24  TYPE(t_field),POINTER :: f_qu(:), f_qv(:), f_Kv(:), f_hv(:), & 
    2525                           f_pk(:),f_wwuu(:),f_planetvel(:), & 
    2626                           f_Fel(:), f_gradPhi2(:), f_wil(:), f_Wetadot(:) 
Note: See TracChangeset for help on using the changeset viewer.