Changeset 377


Ignore:
Timestamp:
04/14/16 23:12:41 (8 years ago)
Author:
dubos
Message:

New : positive advection option for theta

Location:
codes/icosagcm/trunk/src
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed

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

    r373 r377  
    153153       dtheta_rhodz=f_dtheta_rhodz(ind) 
    154154       du=f_du_slow(ind) 
     155 
    155156       IF(hydrostatic) THEN 
    156157          CALL compute_caldyn_slow_hydro(u,mass,hflux,du) 

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

    r373 r377  
    302302       !DIR$ SIMD 
    303303       DO ij=ij_begin,ij_end 
    304           dW(ij,l) = dW(ij,l) + W_etadot(ij,l-1) - W_etadot(ij,l) 
    305304          dPhi(ij,l) = dPhi(ij,l) - wflux(ij,l) & 
    306305               * (geopot(ij,l+1)-geopot(ij,l-1))/(mass(ij,l-1)+mass(ij,l)) 
    307306       END DO 
    308307    END DO 
     308    DO l=ll_begin,ll_end 
     309       !DIR$ SIMD 
     310       DO ij=ij_begin,ij_end 
     311          dW(ij,l+1) = dW(ij,l+1) + W_etadot(ij,l) ! update inner+top interfaces 
     312          dW(ij,l)   = dW(ij,l)   - W_etadot(ij,l) ! update bottom+inner interfaces 
     313       END DO 
     314    END DO 
    309315    CALL trace_end("compute_caldyn_vert_nh") 
    310316 

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

    r375 r377  
    155155 
    156156    ! Newton-Raphson iteration : Phi_il contains current guess value 
    157     DO iter=1,2 ! 2 iterations should be enough 
     157    DO iter=1,5 ! 2 iterations should be enough 
    158158 
    159159       ! Compute pressure, A_ik 
    … …  
    673673          kup=l 
    674674       END IF 
     675       ! below : "checked" means "formula also valid when kup=kdown (top/bottom)" 
    675676       ! compute mil, wil=Wil/mil 
    676677       DO ij=ij_begin_ext, ij_end_ext 
    677           w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup)) 
     678          w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup)) ! checked 
    678679       END DO 
    679680       ! compute DePhi, v_el, G_el, F_el 
    … …  
    687688          W2_el = .5*le_de(ij+u_right) * & 
    688689               ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) ) 
    689           v_el(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) 
     690          v_el(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked 
    690691          G_el(ij+u_right) = v_el(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el 
    691692          ! Compute on edge 'lup' 
    … …  
    695696          W2_el = .5*le_de(ij+u_lup) * & 
    696697               ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) ) 
    697           v_el(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown) ) 
     698          v_el(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked 
    698699          G_el(ij+u_lup) = v_el(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el 
    699700          ! Compute on edge 'ldown' 
    … …  
    703704          W2_el = .5*le_de(ij+u_ldown) * & 
    704705               ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) ) 
    705           v_el(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown) ) 
     706          v_el(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked 
    706707          G_el(ij+u_ldown) = v_el(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el 
    707708       END DO 
    … …  
    715716               le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 +  & 
    716717               le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 ) 
    717           dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l)-1/(2*Ai(ij))* &  
     718!          gradPhi2(ij,l) = 0. ! FIXME !! 
     719 
     720          dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) -1/(2*Ai(ij))* &  
    718721               ( DePhil(ij+u_right,l)*v_el(ij+u_right) + & ! -v.gradPhi,  
    719722                 DePhil(ij+u_rup,l)*v_el(ij+u_rup) +     & ! v_el already has le_de 
    … …  
    722725                 DePhil(ij+u_ldown,l)*v_el(ij+u_ldown) + & 
    723726                 DePhil(ij+u_rdown,l)*v_el(ij+u_rdown) ) 
     727 
    724728          dW(ij,l) = -1./Ai(ij)*(           & ! -div(G_el),  
    725729               ne_right*G_el(ij+u_right) +  & ! G_el already has le_de 
    … …  
    731735       END DO 
    732736    END DO 
    733      
     737    ! FIXME !! 
     738 !   F_el(:,:)=0.  
     739 !   dPhi(:,:)=0. 
     740 !   dW(:,:)=0. 
     741 
    734742    DO l=ll_begin, ll_end ! compute on k levels (layers) 
    735743       ! Compute berni at scalar points 
    … …  
    765773       END DO 
    766774    END DO 
     775    ! FIXME !! 
     776    ! du(:,:)=0. 
     777    ! hflux(:,:)=0. 
    767778 
    768779    CALL trace_end("compute_caldyn_slow_NH") 

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

    r369 r377  
    11961196! Test 3-1 
    11971197!========== 
    1198 SUBROUTINE test3_gravity_wave (lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q) 
     1198SUBROUTINE test3_gravity_wave (X,lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q) 
    11991199 
    12001200IMPLICIT NONE 
    … …  
    12041204 
    12051205        real(rstd), intent(in)  :: lon, &       ! Longitude (radians) 
    1206                                    lat          ! Latitude (radians) 
     1206                                   lat, &       ! Latitude (radians) 
     1207                                   X            ! Reduced Earth reduction factor (DCMIP value = 125) 
    12071208 
    12081209        real(rstd), intent(inout) :: p, &       ! Pressure  (Pa) 
    … …  
    12271228!     test case parameters 
    12281229!-----------------------------------------------------------------------  
    1229         real(rstd), parameter ::        X       = 125.d0,               &       ! Reduced Earth reduction factor 
     1230        real(rstd), parameter ::        & 
    12301231                                Om      = 0.d0,                 &       ! Rotation Rate of Earth 
    1231                                 as      = a/X,                  &       ! New Radius of small Earth      
    12321232                                u0      = 20.d0,                &       ! Reference Velocity  
    12331233!                               u0      = 0.d0,                 &       ! FIXME : no zonal wind for NH tests 
    … …  
    12431243                                N2      = N*N,                  &       ! Brunt-Vaisala frequency Squared 
    12441244                                bigG    = (g*g)/(N2*cp)                 ! Constant 
    1245                              
     1245 
     1246      real(rstd) :: as                                                  ! New Radius of small Earth      
     1247 
    12461248      real(rstd) :: height                                                      ! Model level height 
    12471249      real(rstd) :: sin_tmp, cos_tmp                                    ! Calculation of great circle distance 
    … …  
    12511253      real(rstd) :: theta_pert                                          ! Pot-temp perturbation 
    12521254 
     1255      as = a/X 
     1256 
    12531257!----------------------------------------------------------------------- 
    12541258!    THE VELOCITIES 

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

    r266 r377  
    11MODULE disvert_std_mod 
    22  USE icosa 
     3  IMPLICIT NONE 
     4  PRIVATE 
     5 
    36  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:) 
    47!$OMP THREADPRIVATE(ap) 
    … …  
    710  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:) 
    811!$OMP THREADPRIVATE(presnivs) 
     12 
     13  PUBLIC :: init_disvert, ap, bp, presnivs 
    914 
    1015CONTAINS 

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

    r366 r377  
    1717 
    1818  LOGICAL, SAVE :: boussinesq, hydrostatic 
     19  !$OMP THREADPRIVATE(boussinesq, hydrostatic) 
    1920 
    2021CONTAINS 

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

    r366 r377  
    2121    USE etat0_dcmip4_mod, ONLY : getin_etat0_dcmip4=>getin_etat0 
    2222    USE etat0_dcmip5_mod, ONLY : getin_etat0_dcmip5=>getin_etat0 
     23    USE etat0_bubble_mod, ONLY : getin_etat0_bubble=>getin_etat0 
    2324    USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0 
    2425    USE etat0_temperature_mod, ONLY: getin_etat0_temperature=>getin_etat0 
    … …  
    7172    CASE ('dcmip5') 
    7273        CALL getin_etat0_dcmip5 
     74    CASE ('bubble') 
     75        CALL getin_etat0_bubble 
    7376    CASE ('williamson91.6') 
    7477       autoinit_mass=.FALSE. 
    … …  
    9699      ELSE 
    97100         PRINT*, 'Bad selector for variable etat0 <',etat0_type, & 
    98               '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> ' 
     101              '> options are <isothermal>, <temperature_profile>, <held_suarez>, & 
     102              &<bubble>, <venus>, <jablonowsky06>, <academic>, <dcmip[1-4]> ' 
    99103         STOP 
    100104      END IF 
    … …  
    176180    USE etat0_dcmip4_mod, ONLY : compute_dcmip4 => compute_etat0 
    177181    USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0 
     182    USE etat0_bubble_mod, ONLY : compute_bubble => compute_etat0   
    178183    USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0 
    179184    USE etat0_temperature_mod, ONLY: compute_etat0_temperature => compute_etat0 
    … …  
    205210 
    206211    ! For NH geopotential and vertical momentum must be initialized. 
    207     ! Unless autoinit_mass is set to .FALSE. , they will be initialized 
     212    ! Unless autoinit_NH is set to .FALSE. , they will be initialized 
    208213    ! to hydrostatic geopotential and zero 
    209214    autoinit_mass = .TRUE. 
    … …  
    239244       CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 
    240245       CALL compute_dcmip5(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 
     246    CASE('bubble') 
     247       CALL compute_bubble(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q) 
     248       CALL compute_bubble(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e) 
     249!       autoinit_NH = .FALSE. ! compute_bubble initializes geopot 
    241250    CASE('williamson91.6') 
    242251       CALL compute_w91_6(iim*jjm,lon_i,lat_i, phis, mass(:,1), temp_i(:,1), ulon_i(:,1), ulat_i(:,1)) 

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

    r367 r377  
    88   
    99  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,geopot,q) 
    10     USE genmod 
    1110    USE dcmip_initial_conditions_test_1_2_3 
    1211    USE disvert_mod 
    … …  
    2726    pp=peq 
    2827    DO ij=1,ngrid 
    29        CALL test3_gravity_wave(lon(ij),lat(ij),pp,dummy,0, & 
     28       CALL test3_gravity_wave(scale_factor, lon(ij),lat(ij),pp,dummy,0, & 
    3029            dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy) 
    3130    END DO 
    … …  
    3332       DO ij=1,ngrid 
    3433          pp = ap(l) + bp(l)*ps(ij) ! half-layer pressure 
    35           CALL test3_gravity_wave(lon(ij),lat(ij),pp,zz,0, & 
     34          CALL test3_gravity_wave(scale_factor, lon(ij),lat(ij),pp,zz,0, & 
    3635               dummy,dummy,dummy,dummy,dummy,dummy,dummy,dummy) 
    3736          geopot(ij,l) = g*zz ! initialize geopotential for NH 
    … …  
    4140       DO ij=1,ngrid 
    4241          pp = .5*(ap(l)+ap(l+1)) + .5*(bp(l)+bp(l+1))*ps(ij) ! full-layer pressure 
    43           CALL test3_gravity_wave(lon(ij),lat(ij),pp,dummy,0, & 
     42          CALL test3_gravity_wave(scale_factor, lon(ij),lat(ij),pp,dummy,0, & 
    4443               ulon(ij,l),ulat(ij,l),dummy,Temp(ij,l),dummy,dummy,dummy,dummy) 
    4544       END DO 

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

    r328 r377  
    33  USE icosa, ONLY: llm,nqtot 
    44  IMPLICIT NONE 
     5  PRIVATE 
     6 
    57  REAL(rstd), SAVE, ALLOCATABLE :: t_profile(:) 
    68!$OMP THREADPRIVATE(t_profile) 
     9 
     10  PUBLIC :: getin_etat0, compute_etat0 
     11 
    712CONTAINS 
    813   
    … …  
    1318    USE transfert_omp_mod, ONLY: bcast_omp 
    1419    USE free_unit_mod, ONLY : free_unit 
    15     IMPLICIT NONE 
    1620    INTEGER :: unit,ok 
    1721    INTEGER :: l 
    … …  
    5963  SUBROUTINE compute_etat0(ngrid, phis, ps, temp, ulon, ulat, q) 
    6064    USE earth_const, ONLY: preff 
    61     IMPLICIT NONE   
    6265    INTEGER, INTENT(IN)    :: ngrid 
    6366    REAL(rstd),INTENT(OUT) :: phis(ngrid) 

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

    r374 r377  
    169169    REAL(rstd) :: F_el(3*iim*jjm,llm+1) 
    170170    REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, DePhil 
     171    ! NB : u and uh are not in DEC form, they are normal components     
     172    ! => we must divide by de 
    171173    IF(hydrostatic) THEN 
    172174       uh(:,:)=u(:,:) 
    … …  
    178180             W_el   = .5*( W(ij,l)+W(ij+t_right,l) ) 
    179181             DePhil = ne_right*(Phi(ij+t_right,l)-Phi(ij,l)) 
    180              F_el(ij+u_right,l)   = DePhil*W_el 
     182             F_el(ij+u_right,l)   = DePhil*W_el/de(ij+u_right) 
    181183             ! Compute on edge 'lup' 
    182184             W_el   = .5*( W(ij,l)+W(ij+t_lup,l) ) 
    183185             DePhil = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l)) 
    184              F_el(ij+u_lup,l)   = DePhil*W_el 
     186             F_el(ij+u_lup,l)   = DePhil*W_el/de(ij+u_lup) 
    185187             ! Compute on edge 'ldown' 
    186188             W_el   = .5*( W(ij,l)+W(ij+t_ldown,l) ) 
    187189             DePhil = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l)) 
    188              F_el(ij+u_ldown,l) = DePhil*W_el 
     190             F_el(ij+u_ldown,l) = DePhil*W_el/de(ij+u_ldown) 
    189191          END DO 
    190192       END DO 
    … …  
    192194       DO l=ll_begin, ll_end ! compute on k levels (full levels) 
    193195          DO ij=ij_begin_ext, ij_end_ext 
    194              ! w = vertical momentum = g^-2*dPhi/dt = g*uz 
     196             ! w = vertical momentum = g^-2*dPhi/dt = uz/g 
    195197             uz(ij,l) = (.5*g)*(W(ij,l)+W(ij,l+1))/rhodz(ij,l) 
    196198             ! uh = u-w.grad(Phi) = u - uz.grad(z) 

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

    r376 r377  
    1111   
    1212  INTEGER, PARAMETER :: itau_sync=10 
    13   TYPE(t_message),SAVE :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0 
     13  TYPE(t_message),SAVE :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0, req_W0, req_geopot0 
     14  LOGICAL, SAVE :: positive_theta 
    1415 
    1516  PUBLIC :: init_timeloop, timeloop 
    … …  
    3334    IF (xios_output) itau_out=1 
    3435    IF (.NOT. enable_io) itau_out=HUGE(itau_out) 
    35      
     36 
     37    positive_theta=.FALSE. 
     38    CALL getin('positive_theta',positive_theta) 
     39    IF(positive_theta .AND. nqtot<1) THEN 
     40       PRINT *, 'nqtot must be >0 if positive_theta is .TRUE.' 
     41       STOP 
     42    END IF 
     43 
    3644    def='ARK2.3' 
    3745    CALL getin('time_scheme',def) 
    … …  
    147155    CALL init_message(f_u,req_e0_vect,req_u0) 
    148156    CALL init_message(f_q,req_i0,req_q0) 
     157    CALL init_message(f_geopot,req_i0,req_geopot0) 
     158    CALL init_message(f_W,req_i0,req_W0) 
    149159 
    150160  END SUBROUTINE init_timeloop 
    … …  
    198208    fluxt_zero=.TRUE. 
    199209 
     210    IF(positive_theta) CALL copy_theta_to_q1(f_theta_rhodz,f_rhodz,f_q) 
     211 
    200212    !$OMP MASTER 
    201213    CALL SYSTEM_CLOCK(start_clock) 
    … …  
    226238          CALL wait_message(req_u0) 
    227239          CALL send_message(f_q,req_q0)  
    228           CALL wait_message(req_q0)  
     240          CALL wait_message(req_q0) 
     241          IF(.NOT. hydrostatic) THEN 
     242!             CALL send_message(f_geopot,req_geopot0) 
     243!             CALL wait_message(req_geopot0) 
     244!             CALL send_message(f_W,req_W0) 
     245!             CALL wait_message(req_W0) 
     246          END IF 
    229247       ENDIF 
    230248 
    … …  
    285303             END DO 
    286304          ENDIF 
     305          IF(positive_theta) CALL copy_q1_to_theta(f_theta_rhodz,f_rhodz,f_q) 
    287306       END IF 
    288307        
    … …  
    341360  END SUBROUTINE print_iteration 
    342361 
     362  SUBROUTINE copy_theta_to_q1(f_theta_rhodz,f_rhodz,f_q) 
     363    TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:) 
     364    REAL(rstd), POINTER :: theta_rhodz(:,:), rhodz(:,:), q(:,:,:) 
     365    INTEGER :: ind 
     366    DO ind=1, ndomain 
     367       IF (.NOT. assigned_domain(ind)) CYCLE 
     368       CALL swap_dimensions(ind) 
     369       CALL swap_geometry(ind) 
     370       theta_rhodz=f_theta_rhodz(ind) 
     371       rhodz=f_rhodz(ind) 
     372       q=f_q(ind) 
     373       q(:,:,1)  = theta_rhodz(:,:)/rhodz(:,:) 
     374    END DO 
     375  END SUBROUTINE copy_theta_to_q1 
     376 
     377  SUBROUTINE copy_q1_to_theta(f_theta_rhodz,f_rhodz,f_q) 
     378    TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:) 
     379    REAL(rstd), POINTER :: theta_rhodz(:,:), rhodz(:,:), q(:,:,:) 
     380    INTEGER :: ind 
     381    DO ind=1, ndomain 
     382       IF (.NOT. assigned_domain(ind)) CYCLE 
     383       CALL swap_dimensions(ind) 
     384       CALL swap_geometry(ind) 
     385       theta_rhodz=f_theta_rhodz(ind) 
     386       rhodz=f_rhodz(ind) 
     387       q=f_q(ind) 
     388       theta_rhodz(:,:) = rhodz(:,:)*q(:,:,1) 
     389    END DO 
     390  END SUBROUTINE copy_q1_to_theta 
     391 
    343392END MODULE timeloop_gcm_mod 
Note: See TracChangeset for help on using the changeset viewer.