Changeset 366 for codes


Ignore:
Timestamp:
10/30/15 15:41:06 (9 years ago)
Author:
dubos
Message:

Progress towards NH

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

Legend:

Unmodified
Added
Removed

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

    r362 r366  
    1212   
    1313  SUBROUTINE caldyn_hevi(write_out,tau, f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & 
    14        f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du_slow, f_du_fast) 
     14       f_W, f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, & 
     15       f_du_slow, f_du_fast, f_dPhi_slow, f_dPhi_fast, f_dW_slow, f_dW_fast)  
    1516    USE icosa 
    1617    USE observable_mod 
    … …  
    3435    TYPE(t_field),POINTER :: f_u(:) 
    3536    TYPE(t_field),POINTER :: f_q(:) 
     37    TYPE(t_field),POINTER :: f_W(:) 
    3638    TYPE(t_field),POINTER :: f_geopot(:) 
    3739    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:) 
    … …  
    4143    TYPE(t_field) :: f_du_slow(:) 
    4244    TYPE(t_field) :: f_du_fast(:) 
     45    TYPE(t_field) :: f_dW_slow(:) 
     46    TYPE(t_field) :: f_dW_fast(:) 
     47    TYPE(t_field) :: f_dPhi_slow(:) 
     48    TYPE(t_field) :: f_dPhi_fast(:) 
    4349     
    4450    REAL(rstd),POINTER :: ps(:), dps(:) 
    4551    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:) 
    46     REAL(rstd),POINTER :: du(:,:), hflux(:,:), wflux(:,:) 
    47     REAL(rstd),POINTER :: u(:,:), qu(:,:), qv(:,:) 
     52    REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:) 
     53    REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:) 
    4854 
    4955! temporary shared variable 
    … …  
    6874      CALL init_message(f_u,req_e1_vect,req_u) 
    6975      CALL init_message(f_qu,req_e1_scal,req_qu) 
     76      IF(.NOT.hydrostatic) THEN 
     77         CALL init_message(f_geopot,req_i1,req_geopot) 
     78         CALL init_message(f_w,req_i1,req_w) 
     79      END IF 
    7080    ENDIF 
    7181     
    … …  
    8191    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01 
    8292    CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort 
     93 
     94    IF(.NOT.hydrostatic) THEN 
     95       CALL send_message(f_geopot,req_geopot) ! COM03 
     96       CALL wait_message(req_geopot) ! COM03 
     97       CALL send_message(f_w,req_w) ! COM04 
     98       CALL wait_message(req_w) ! COM04 
     99    END IF 
    83100     
    84101    DO ind=1,ndomain 
    … …  
    87104       CALL swap_geometry(ind) 
    88105       ps=f_ps(ind) 
    89        u=f_u(ind) 
    90        du=f_du_fast(ind) 
    91106       theta_rhodz=f_theta_rhodz(ind) 
    92107       mass=f_mass(ind) 
    93108       theta = f_theta(ind) 
     109       CALL compute_theta(ps,theta_rhodz, mass,theta) 
    94110       pk = f_pk(ind) 
    95111       geopot = f_geopot(ind) 
    96        CALL compute_theta(ps,theta_rhodz, mass,theta) 
    97        CALL compute_geopot(ps,mass,theta, pk,geopot) 
    98        CALL compute_caldyn_fast(tau,u,mass,theta,pk,geopot, du) ! computes du_fast and updates u 
     112       IF(hydrostatic) THEN 
     113          CALL compute_geopot(ps,mass,theta, pk,geopot) 
     114       ELSE 
     115          W = f_W(ind) 
     116          dW = f_dW_fast(ind) 
     117          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 
     119       END IF 
     120       u=f_u(ind) 
     121       du=f_du_fast(ind) 
     122       CALL compute_caldyn_fast(tau,u,mass,theta,pk,geopot,du) ! computes du_fast and updates du 
    99123    ENDDO 
    100124     

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

    r362 r366  
    1818  !$OMP THREADPRIVATE(caldyn_conserv)  
    1919 
    20   TYPE(t_message),PUBLIC :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu 
     20  TYPE(t_message),PUBLIC :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu, req_geopot, req_w 
    2121 
    2222  PUBLIC :: compute_geopot, compute_caldyn_vert 
    … …  
    5151 
    5252    IF(caldyn_eta==eta_mass .AND. .NOT. DEC) THEN 
    53  
    5453!!! Compute exner function and geopotential 
    5554       DO l = 1,llm 
    … …  
    6463       ENDDO 
    6564    ELSE  
    66        ! We are using a Lagrangian vertical coordinate 
    67        ! Pressure must be computed first top-down (temporarily stored in pk) 
     65       ! We are using DEC or a Lagrangian vertical coordinate 
     66       ! Pressure is computed first top-down (temporarily stored in pk) 
    6867       ! Then Exner pressure and geopotential are computed bottom-up 
    69        ! Notice that the computation below should work also when caldyn_eta=eta_mass           
     68       ! Works also when caldyn_eta=eta_mass           
    7069 
    7170       IF(boussinesq) THEN ! compute geopotential and pk=Lagrange multiplier 

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

    r362 r366  
    66  PRIVATE 
    77 
    8   PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_fast, compute_caldyn_slow 
     8  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_slow, & 
     9       compute_caldyn_solver, compute_caldyn_fast 
    910 
    1011CONTAINS 
    … …  
    120121  END SUBROUTINE compute_pvort_only 
    121122 
    122   SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot, du) 
     123  SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk, geopot,W, dPhi,dW) 
     124    USE icosa 
     125    USE disvert_mod 
     126    USE exner_mod 
     127    USE trace 
     128    USE omp_para 
     129    USE disvert_mod 
     130    IMPLICIT NONE 
     131    REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs 
     132    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm) 
     133    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm) 
     134    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm) 
     135    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) 
     136    REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0 
     137    REAL(rstd),INTENT(OUT)   :: dW(iim*jjm,llm+1) 
     138    REAL(rstd),INTENT(OUT)   :: dPhi(iim*jjm,llm+1) 
     139 
     140    INTEGER :: ij,l 
     141    INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext 
     142    REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij, m_il 
     143 
     144    CALL trace_start("compute_caldyn_solver") 
     145 
     146    CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext) 
     147    ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1 
     148    ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1 
     149 
     150    IF(tau>0) THEN 
     151       PRINT *, 'No implicit solver for NH yet' 
     152       STOP 
     153    END IF 
     154 
     155    ! Compute pressure, stored temporarily in pk 
     156    ! kappa = R/Cp 
     157    ! 1-kappa = Cv/Cp 
     158    ! Cp/Cv = 1/(1-kappa) 
     159    gamma = 1./(1.-kappa) 
     160    DO l=ll_begin,ll_end 
     161       !DIR$ SIMD 
     162       DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     163          rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l)) 
     164          X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij 
     165          ! kappa.theta.rho = p/exner  
     166          ! => X = (p/p0)/(exner/Cp) 
     167          !      = (p/p0)^(1-kappa) 
     168          pk(ij,l) = preff*(X_ij**gamma) 
     169       ENDDO 
     170    ENDDO 
     171 
     172    ! Compute tendencies for geopotential and vertical momentum 
     173    DO l=ll_beginp1,ll_end 
     174       !DIR$ SIMD 
     175       DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     176          m_il = .5*(rhodz(ij,l)+rhodz(ij,l-1)) 
     177          dW(ij,l) = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il 
     178          dPhi(ij,l) = g*g*W(ij,l)/m_il 
     179       ENDDO 
     180       !       PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l))) 
     181       !       PRINT *,'Max dW', l,ij_begin,ij_end, MAXVAL(abs(dW(ij_begin:ij_end,l)))  
     182    ENDDO 
     183    ! Lower BC (FIXME : no orography yet !) 
     184    DO ij=ij_begin,ij_end          
     185       dPhi(ij,1)=0 
     186       W(ij,1)=0 
     187       dW(ij,1)=0 
     188       dPhi(ij,llm+1)=0 ! rigid lid 
     189       W(ij,llm+1)=0 
     190       dW(ij,llm+1)=0 
     191    ENDDO 
     192    ! Upper BC p=ptop 
     193!    DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     194!       dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm) 
     195!       dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm) 
     196!    ENDDO 
     197     
     198    ! Compute Exner function (needed by compute_caldyn_fast) 
     199    DO l=ll_begin,ll_end 
     200       !DIR$ SIMD 
     201       DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     202          pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow 
     203       ENDDO 
     204    ENDDO 
     205     
     206    CALL trace_end("compute_caldyn_solver") 
     207     
     208  END SUBROUTINE compute_caldyn_solver 
     209   
     210  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du) 
    123211    USE icosa 
    124212    USE disvert_mod 
    … …  
    127215    USE omp_para 
    128216    IMPLICIT NONE 
    129     REAL(rstd), INTENT(IN) :: tau ! "solve" u-tau*du/dt = rhs 
    130     REAL(rstd),INTENT(INOUT)  :: u(iim*3*jjm,llm)    ! prognostic "velocity" 
    131     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    132     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature 
    133     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function 
    134     REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential 
    135     REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) 
     217    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs 
     218    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! OUT if tau>0 
     219    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm) 
     220    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm) 
     221    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) 
     222    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) 
     223    REAL(rstd),INTENT(OUT)   :: du(iim*3*jjm,llm) 
    136224    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
    137225 
    … …  
    140228 
    141229    CALL trace_start("compute_caldyn_fast") 
    142      
    143     ! Compute bernouilli term 
     230 
     231    ! Compute Bernoulli term 
    144232    IF(boussinesq) THEN 
     233 
    145234       DO l=ll_begin,ll_end 
    146235          !DIR$ SIMD 

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

    r186 r366  
    12031203!----------------------------------------------------------------------- 
    12041204 
    1205         real(rstd), intent(in)  :: lon, &               ! Longitude (radians) 
    1206                                 lat, &          ! Latitude (radians) 
    1207                                 z               ! Height (m) 
    1208  
    1209         real(rstd), intent(inout) :: p          ! Pressure  (Pa) 
    1210                                  
     1205        real(rstd), intent(in)  :: lon, &       ! Longitude (radians) 
     1206                                   lat          ! Latitude (radians) 
     1207 
     1208        real(rstd), intent(inout) :: p, &       ! Pressure  (Pa) 
     1209                                     z          ! Height (m) 
     1210 
    12111211 
    12121212        integer,  intent(in) :: zcoords         ! 0 or 1 see below 
    12131213 
    1214         real(rstd), intent(out) :: u, &                 ! Zonal wind (m s^-1) 
     1214        real(rstd), intent(out) :: u, &         ! Zonal wind (m s^-1) 
    12151215                                v, &            ! Meridional wind (m s^-1) 
    12161216                                w, &            ! Vertical Velocity (m s^-1) 
    … …  
    12301230                                Om      = 0.d0,                 &       ! Rotation Rate of Earth 
    12311231                                as      = a/X,                  &       ! New Radius of small Earth      
    1232                                 u0      = 20.d0,                &       ! Reference Velocity  
     1232!                               u0      = 20.d0,                &       ! Reference Velocity  
     1233                                u0      = 0.d0,                 &       ! FIXME : no zonal wind for NH tests 
    12331234                                Teq     = 300.d0,               &       ! Temperature at Equator     
    12341235                                Peq     = 100000.d0,            &       ! Reference PS at Equator 
    … …  
    12971298 
    12981299                height = (-g/N2)*log( (TS/bigG)*( (p/ps)**(Rd/cp) - 1.d0  ) + 1.d0 ) 
     1300                z = height ! modified : initialize z when pressure-based 
    12991301 
    13001302        endif 

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

    r266 r366  
    1616  REAL(rstd),SAVE :: mu                 ! molar mass of the atmosphere 
    1717 
    18   LOGICAL, SAVE :: boussinesq 
     18  LOGICAL, SAVE :: boussinesq, hydrostatic 
    1919 
    2020CONTAINS 
    … …  
    3434    CALL getin("scale_height",scale_height) 
    3535     
    36     mu=kappa/cpp 
    3736    boussinesq=.FALSE. 
    3837    CALL getin("boussinesq",boussinesq)   
     38    hydrostatic=.TRUE. 
     39    CALL getin("hydrostatic",hydrostatic)   
     40    IF(boussinesq .AND. .NOT. hydrostatic) THEN 
     41       PRINT *, 'boussinesq=.TRUE. and hydrostatic=.FALSE. : Non-hydrostatic boussinesq equations are not supported' 
     42       STOP 
     43    END IF 
    3944     
     45    mu=kappa/cpp 
    4046    radius=radius/scale_factor 
    4147    omega=omega*scale_factor 
    … …  
    4551   
    4652END MODULE earth_const 
    47    

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

    r353 r366  
    1313CONTAINS 
    1414   
    15   SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q) 
     15  SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_w, f_q) 
    1616    USE mpipara, ONLY : is_mpi_root 
    1717    USE disvert_mod 
    … …  
    3535    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    3636    TYPE(t_field),POINTER :: f_u(:) 
     37    TYPE(t_field),POINTER :: f_geopot(:) 
     38    TYPE(t_field),POINTER :: f_w(:) 
    3739    TYPE(t_field),POINTER :: f_q(:) 
    3840     
    3941    REAL(rstd),POINTER :: ps(:), mass(:,:) 
    40     LOGICAL :: init_mass, collocated 
     42    LOGICAL :: autoinit_mass, autoinit_geopot, collocated 
    4143    INTEGER :: ind,i,j,ij,l 
    4244 
    … …  
    4547    ! the initial distribution of mass is taken to be the same 
    4648    ! as what the mass coordinate would dictate 
    47     ! however if etat0_XXX defines mass then the flag init_mass must be set to .FALSE. 
     49    ! however if etat0_XXX defines mass then the flag autoinit_mass must be set to .FALSE. 
    4850    ! otherwise mass will be overwritten  
    49     init_mass = (caldyn_eta == eta_lag) 
     51    autoinit_mass = (caldyn_eta == eta_lag) 
    5052 
    5153    etat0_type='jablonowsky06' 
    … …  
    7072        CALL getin_etat0_dcmip5 
    7173    CASE ('williamson91.6') 
    72        init_mass=.FALSE. 
     74       autoinit_mass=.FALSE. 
    7375       CALL getin_etat0_williamson 
    7476    CASE DEFAULT 
    7577       collocated=.FALSE. 
     78       autoinit_mass = .FALSE. 
    7679    END SELECT 
    7780 
    … …  
    9093   CASE DEFAULT 
    9194      IF(collocated) THEN 
    92          CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q) 
     95         CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q) 
    9396      ELSE 
    9497         PRINT*, 'Bad selector for variable etat0 <',etat0_type, & 
    … …  
    98101    END SELECT 
    99102 
    100     IF(init_mass) THEN ! initialize mass distribution using ps 
    101103!       !$OMP BARRIER 
     104    IF(autoinit_mass) THEN 
    102105       DO ind=1,ndomain 
    103106          IF (.NOT. assigned_domain(ind)) CYCLE 
    … …  
    105108          CALL swap_geometry(ind) 
    106109          mass=f_mass(ind); ps=f_ps(ind) 
    107           CALL compute_rhodz(.TRUE., ps, mass) 
     110          CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps 
    108111       END DO 
    109112    END IF 
    110  
     113  
    111114  END SUBROUTINE etat0 
    112115 
    113   SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q) 
     116  SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q) 
    114117    USE theta2theta_rhodz_mod 
    115118    IMPLICIT NONE 
    … …  
    119122    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    120123    TYPE(t_field),POINTER :: f_u(:) 
     124    TYPE(t_field),POINTER :: f_geopot(:) 
     125    TYPE(t_field),POINTER :: f_W(:) 
    121126    TYPE(t_field),POINTER :: f_q(:) 
    122127   
    … …  
    128133    REAL(rstd),POINTER :: temp(:,:) 
    129134    REAL(rstd),POINTER :: u(:,:) 
     135    REAL(rstd),POINTER :: geopot(:,:) 
     136    REAL(rstd),POINTER :: W(:,:) 
    130137    REAL(rstd),POINTER :: q(:,:,:) 
    131138    INTEGER :: ind 
    … …  
    143150      temp=f_temp(ind) 
    144151      u=f_u(ind) 
     152      geopot=f_geopot(ind) 
     153      w=f_w(ind) 
    145154      q=f_q(ind) 
    146155 
    147       IF( TRIM(etat0_type)=='williamson91.6' ) THEN  
    148        CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q) 
     156      IF( TRIM(etat0_type)=='williamson91.6' ) THEN 
     157       CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, geopot, W, q) 
    149158      ELSE 
    150        CALL compute_etat0_collocated(ps,mass, phis, temp, u, q) 
     159       CALL compute_etat0_collocated(ps,mass, phis, temp, u, geopot, W, q) 
    151160      ENDIF 
    152161    ENDDO 
    … …  
    158167  END SUBROUTINE etat0_collocated 
    159168 
    160   SUBROUTINE compute_etat0_collocated(ps,mass, phis, temp_i, u, q) 
     169  SUBROUTINE compute_etat0_collocated(ps,mass,phis,temp_i,u, geopot,W, q) 
    161170    USE wind_mod 
     171    USE disvert_mod 
    162172    USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0 
    163173    USE etat0_dcmip1_mod, ONLY : compute_dcmip1 => compute_etat0 
    … …  
    174184    REAL(rstd),INTENT(OUT) :: temp_i(iim*jjm,llm) 
    175185    REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm) 
     186    REAL(rstd),INTENT(OUT) :: W(iim*jjm,llm+1) 
     187    REAL(rstd),INTENT(OUT) :: geopot(iim*jjm,llm+1) 
    176188    REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot) 
    177189 
    … …  
    183195    REAL(rstd) :: phis_e(3*iim*jjm) 
    184196    REAL(rstd) :: temp_e(3*iim*jjm,llm) 
     197    REAL(rstd) :: geopot_e(3*iim*jjm,llm+1) 
    185198    REAL(rstd) :: ulon_e(3*iim*jjm,llm) 
    186199    REAL(rstd) :: ulat_e(3*iim*jjm,llm) 
    … …  
    188201 
    189202    INTEGER :: l,i,j,ij 
     203    REAL :: p_ik, v_ik, mass_ik 
     204    LOGICAL :: autoinit_mass, autoinit_NH 
     205 
     206    ! For NH geopotential and vertical momentum must be initialized. 
     207    ! Unless autoinit_mass is set to .FALSE. , they will be initialized 
     208    ! to hydrostatic geopotential and zero 
     209    autoinit_mass = .TRUE. 
     210    autoinit_NH = .NOT. hydrostatic 
     211    w(:,:) = 0 
    190212 
    191213    !$OMP BARRIER 
    … …  
    208230       CALL compute_dcmip2(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)       
    209231    CASE('dcmip3') 
    210        CALL compute_dcmip3(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 
    211        CALL compute_dcmip3(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 
     232       CALL compute_dcmip3(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q) 
     233       CALL compute_dcmip3(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e) 
     234       autoinit_NH = .FALSE. ! compute_dcmip3 initializes geopot 
    212235    CASE('dcmip4') 
    213236       CALL compute_dcmip4(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 
    … …  
    219242       CALL compute_w91_6(iim*jjm,lon_i,lat_i, phis, mass(:,1), temp_i(:,1), ulon_i(:,1), ulat_i(:,1)) 
    220243       CALL compute_w91_6(3*iim*jjm,lon_e,lat_e, phis_e, mass_e(:,1), temp_e(:,1), ulon_e(:,1), ulat_e(:,1)) 
     244       autoinit_mass = .FALSE. ! do not overwrite mass 
    221245    END SELECT 
     246 
     247    IF(autoinit_mass) CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps 
     248    IF(autoinit_NH) THEN 
     249       geopot(:,1) = phis(:) ! surface geopotential 
     250       DO l = 1, llm 
     251          DO ij=1,iim*jjm 
     252             ! hybrid pressure coordinate 
     253             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) 
     254             mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g 
     255             ! v=R.T/p, R=kappa*cpp 
     256             v_ik = kappa*cpp*temp_i(ij,l)/p_ik 
     257             geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g 
     258          END DO 
     259       END DO 
     260    END IF 
    222261 
    223262    !$OMP BARRIER 

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

    r356 r366  
    77CONTAINS 
    88   
    9   SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,q) 
     9  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,geopot,q) 
    1010    USE genmod 
    1111    USE dcmip_initial_conditions_test_1_2_3 
    … …  
    2020    REAL(rstd), INTENT(OUT) :: ulat(ngrid,llm) 
    2121    REAL(rstd), INTENT(OUT) :: temp(ngrid,llm) 
     22    REAL(rstd), INTENT(OUT) :: geopot(ngrid,llm+1) 
    2223    REAL(rstd), INTENT(OUT) :: q(ngrid,llm,nqtot) 
    2324    REAL(rstd),PARAMETER :: Peq=1e5        ! Reference surface pressure at the equator (hPa) 
    24     REAL(rstd) :: dummy, pp 
     25    REAL(rstd) :: dummy, pp, zz 
    2526    INTEGER :: l,ij 
    2627    pp=peq 
    … …  
    2829       CALL test3_gravity_wave(lon(ij),lat(ij),pp,dummy,0, & 
    2930            dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy) 
     31    END DO 
     32    DO l=ll_begin,ll_endp1 
     33       DO ij=1,ngrid 
     34          pp = ap(l) + bp(l)*ps(ij) ! half-layer pressure 
     35          CALL test3_gravity_wave(lon(ij),lat(ij),pp,zz,0, & 
     36               dummy,dummy,dummy,dummy,dummy,dummy,dummy,dummy) 
     37          geopot(ij,l) = g*zz ! initialize geopotential for NH 
     38       END DO 
     39       q(:,l,:)=0. 
    3040    END DO 
    3141    DO l=ll_begin,ll_end 

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

    r362 r366  
    55 
    66  TYPE(t_field),POINTER,SAVE,PUBLIC :: f_geopot(:), f_q(:), & 
    7        f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:), & 
    8        f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:), & 
    9        f_u(:),f_um1(:),f_um2(:), f_du(:), & ! previous state 
    10        f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:), & ! trends 
    11        f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:), & ! accumulated mass fluxes 
    12                                 ! slow/fast trend arrays for HEVI time scheme 
    13        f_dmass_slow(:,:), f_dps_slow(:,:), f_dtheta_rhodz_slow(:,:), & 
    14        f_du_slow(:,:), f_du_fast(:,:) 
     7       f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:), & ! current and previous time steps + tendency of mass, 
     8       f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:), &           ! column-integrated mass 
     9       f_u(:),f_um1(:),f_um2(:), f_du(:), f_W(:), &                   ! 'horizontal' and vertical (NH) momentum 
     10       f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:), & ! mass-weighted theta/entropy 
     11       f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:), &            ! accumulated mass fluxes 
     12       f_du_slow(:,:), f_du_fast(:,:), &                              ! slow/fast trend arrays  
     13       f_dmass_slow(:,:), f_dps_slow(:,:), f_dtheta_rhodz_slow(:,:), &! for HEVI time scheme  
     14       f_dPhi_slow(:,:), f_dPhi_fast(:,:), &                          ! geopotential tendencies (NH) 
     15       f_dW_slow(:,:), f_dW_fast(:,:)                                 ! vertical momentum tendencies (NH) 
    1516 
    1617  INTEGER, PARAMETER, PUBLIC :: explicit=1, hevi=2, euler=1, rk4=2, mlf=3, rk25=4, ark23=6, ark33=7 

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

    r362 r366  
    2828  SUBROUTINE set_coefs_ark33(dt) 
    2929    ! Fully-explicit RK3 scheme disguised as ARK 
    30     ! To check correctness of caldyn_hevi 
    3130    REAL(rstd) :: dt 
    32     REAL(rstd), PARAMETER, DIMENSION(3,4) :: & 
    33          ajl=(/ zero, (/.5,0.,0./), (/-1.,2.,0./), (/1./6.,2./3.,1./6./) /) 
    34     CALL set_coefs_hevi(dt, ajl, ajl) 
     31    CALL set_coefs_rk(dt, (/ zero, (/.5,0.,0./), (/-1.,2.,0./), (/1./6.,2./3.,1./6./) /) ) 
    3532  END SUBROUTINE set_coefs_ark33 
    3633     
     34  SUBROUTINE set_coefs_rk(dt, ajl) 
     35    REAL(rstd) :: dt, ajl(3,4) 
     36    CALL set_coefs_hevi(dt,ajl,ajl) 
     37  END SUBROUTINE set_coefs_rk 
     38 
    3739  SUBROUTINE set_coefs_hevi(dt, ajl_slow, ajl_fast) 
    3840    REAL(rstd) :: dt, ajl_slow(3,4), ajl_fast(3,4) ! fast/slow Butcher tableaus 
    … …  
    5860       CALL caldyn_hevi((j==1) .AND. (MOD(it,itau_out)==0), taujj(j), & 
    5961            f_phis, f_ps,f_mass,f_theta_rhodz,f_u,f_q, & 
    60             f_geopot, f_hflux, f_wflux, & 
     62            f_W, f_geopot, f_hflux, f_wflux, & 
    6163            f_dps_slow(:,j), f_dmass_slow(:,j), f_dtheta_rhodz_slow(:,j), & 
    62             f_du_slow(:,j), f_du_fast(:,j) ) 
     64            f_du_slow(:,j), f_du_fast(:,j), & 
     65            f_dPhi_slow(:,j), f_dPhi_fast(:,j), & 
     66            f_dW_slow(:,j), f_dW_fast(:,j) ) 
    6367       ! accumulate mass fluxes for transport scheme 
    6468       DO ind=1,ndomain 
    … …  
    7983          CALL update(bjl(l,j), f_u, f_du_slow(:,l)) 
    8084          CALL update(cjl(l,j), f_u, f_du_fast(:,l)) 
     85          IF(.NOT. hydrostatic) THEN 
     86 !            CALL update(bjl(l,j), f_W, f_dW_slow(:,l)) ! slow tendencies of w, Phi not implemented yet 
     87             CALL update(cjl(l,j), f_W, f_dW_fast(:,l)) 
     88 !            CALL update(bjl(l,j), f_geopot, f_dPhi_slow(:,l)) 
     89             CALL update(cjl(l,j), f_geopot, f_dPhi_fast(:,l)) 
     90          END IF 
    8191       END DO 
    8292    END DO 

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

    r365 r366  
    7676    CALL allocate_field(f_phis,field_t,type_real,name='phis') 
    7777    ! Model state at current time step 
    78     CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot') 
    7978    CALL allocate_field(f_ps,field_t,type_real, name='ps') 
    8079    CALL allocate_field(f_mass,field_t,type_real,llm,name='mass') 
    … …  
    8281    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,name='theta_rhodz') 
    8382    CALL allocate_field(f_u,field_u,type_real,llm,name='u') 
     83    CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot') 
     84    CALL allocate_field(f_W,field_t,type_real,llm+1,name='W') 
    8485    CALL allocate_field(f_q,field_t,type_real,llm,nqtot,'q') 
    8586    ! Mass fluxes 
    … …  
    108109       CALL allocate_fields(nb_stage,f_du_slow, field_u,type_real,llm,name='du_slow') 
    109110       CALL allocate_fields(nb_stage,f_du_fast, field_u,type_real,llm,name='du_fast') 
     111       CALL allocate_fields(nb_stage,f_dW_slow, field_t,type_real,llm+1,name='dW_slow') 
     112       CALL allocate_fields(nb_stage,f_dW_fast, field_t,type_real,llm+1,name='dW_fast') 
     113       CALL allocate_fields(nb_stage,f_dPhi_slow, field_t,type_real,llm+1,name='dPhi_slow') 
     114       CALL allocate_fields(nb_stage,f_dPhi_fast, field_t,type_real,llm+1,name='dPhi_fast') 
    110115       f_dps => f_dps_slow(:,1) 
    111116       f_du => f_du_slow(:,1) 
    … …  
    131136    CALL init_check_conserve 
    132137 
    133     CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q) 
     138    CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_W, f_q) 
    134139 
    135140    CALL transfert_request(f_phis,req_i0)  
Note: See TracChangeset for help on using the changeset viewer.