Changeset 741


Ignore:
Timestamp:
09/28/18 12:59:46 (6 years ago)
Author:
dubos
Message:

devel : added vertical diffusion to idealized venus physics

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

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/devel/src/initial/etat0_venus.f90

    r531 r741  
    33  IMPLICIT NONE 
    44  PRIVATE 
     5  SAVE 
    56 
    67  TYPE(t_field),POINTER :: f_temp_eq( :) 
    78  TYPE(t_field),POINTER :: f_temp(:) ! buffer used for physics 
    8  
    9   REAL(rstd), SAVE :: kfrict 
    10 !$OMP THREADPRIVATE(kfrict) 
    11  
    12   PUBLIC :: etat0, init_physics, physics 
     9  REAL(rstd), ALLOCATABLE :: temp_eq_packed(:,:) 
     10 
     11  REAL(rstd) :: kfrict, kv 
     12!$OMP THREADPRIVATE(kfrict, kv) 
     13 
     14  REAL(rstd), PARAMETER :: tauCLee=86400*25 ! 25 Earth days, cf Lebonnois 2012 
     15 
     16  PUBLIC :: etat0, init_physics, full_physics 
    1317   
    1418CONTAINS 
    … …  
    1620!-------------------------------- "Physics" ---------------------------------------- 
    1721 
    18   SUBROUTINE init_physics 
     22  SUBROUTINE init_physics_old 
    1923    USE getin_mod 
    20     IMPLICIT NONE 
    2124    REAL(rstd),POINTER :: temp(:,:) 
    2225    REAL(rstd) :: friction_time 
    2326    INTEGER :: ind 
    2427 
     28    kv=0.15  ! vertical turbulent viscosity 
     29    CALL getin('venus_diffusion',kv) 
    2530    friction_time=86400.  !friction  
    26     CALL getin('friction_time',friction_time) 
     31    CALL getin('venus_friction_time',friction_time) 
    2732    kfrict=1./friction_time 
    2833 
    29     CALL allocate_field(f_temp,field_t,type_real,llm) ! Buffer for later use 
     34    CALL allocate_field(f_temp,field_t,type_real,llm) ! Buffer for later use by physics 
    3035 
    3136    PRINT *, 'Initializing Temp_eq (venus)' 
    … …  
    3641       CALL swap_geometry(ind) 
    3742       temp=f_temp_eq(ind) 
    38        CALL compute_temp_ref(temp, .TRUE.) ! FIXME With meridional gradient 
    39     ENDDO 
    40   END SUBROUTINE init_physics 
     43       CALL compute_temp_ref(.TRUE.,iim*jjm,lat_i, temp) ! With meridional gradient 
     44    ENDDO 
     45  END SUBROUTINE init_physics_old 
    4146 
    4247  SUBROUTINE physics(f_ps,f_theta_rhodz,f_u)  
    43     USE icosa 
    4448    USE theta2theta_rhodz_mod 
    45     IMPLICIT NONE 
    4649    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    4750    TYPE(t_field),POINTER :: f_u(:) 
    … …  
    6669 
    6770  SUBROUTINE compute_physics(temp_eq, temp, u) 
    68     USE icosa 
    69     USE theta2theta_rhodz_mod 
    70     IMPLICIT NONE 
    7171    REAL(rstd),INTENT(IN)    :: temp_eq(iim*jjm,llm)  
    7272    REAL(rstd),INTENT(INOUT) :: temp(iim*jjm,llm) 
    7373    REAL(rstd),INTENT(INOUT) :: u(3*iim*jjm,llm) 
    74     REAL(rstd), PARAMETER :: tauCLee=86400*25 ! 25 Earth days, cf Lebonnois 2012 
    7574    INTEGER :: i,j,l,ij 
    7675    DO l=1,llm 
    … …  
    8584  END SUBROUTINE compute_physics 
    8685 
     86!----------------------- Re-implementation using physics_interface_mod ----------------- 
     87 
     88  SUBROUTINE init_physics 
     89    USE getin_mod 
     90    USE physics_interface_mod 
     91    REAL(rstd),POINTER :: temp(:,:) 
     92    REAL(rstd) :: friction_time 
     93    INTEGER :: ngrid 
     94 
     95    kv=0.15  ! vertical turbulent viscosity 
     96    CALL getin('venus_diffusion',kv) 
     97    friction_time=86400.  !friction  
     98    CALL getin('venus_friction_time',friction_time) 
     99    kfrict=1./friction_time 
     100 
     101    ngrid = physics_inout%ngrid 
     102    ALLOCATE(temp_eq_packed(ngrid,llm)) 
     103    CALL compute_temp_ref(.TRUE.,ngrid,physics_inout%lat, temp_eq_packed) ! With meridional gradient 
     104  END SUBROUTINE init_physics 
     105 
     106  SUBROUTINE full_physics 
     107    USE physics_interface_mod 
     108    CALL compute_physics_column(physics_inout%ngrid, physics_inout%dt_phys, & 
     109         physics_inout%p, physics_inout%geopot, & 
     110         physics_inout%Temp, physics_inout%ulon, physics_inout%ulat, & 
     111         physics_inout%dTemp, physics_inout%dulon, physics_inout%dulat) 
     112  END SUBROUTINE full_physics 
     113 
     114  SUBROUTINE compute_physics_column(ngrid,dt_phys,p,geopot,Temp,u,v,dTemp,du,dv) 
     115    USE earth_const, only : g 
     116    INTEGER, INTENT(IN)   :: ngrid 
     117    REAL(rstd),INTENT(IN) :: dt_phys, & 
     118         p(ngrid,llm+1), geopot(ngrid,llm+1), & 
     119         Temp(ngrid,llm), u(ngrid,llm), v(ngrid,llm) 
     120    REAL(rstd),INTENT(OUT) :: dTemp(ngrid,llm), du(ngrid,llm), dv(ngrid,llm) 
     121    REAL(rstd) :: mass(ngrid,llm), &  ! rho.dz 
     122                  A(ngrid,llm+1),  &  ! off-diagonal coefficients 
     123                  B(ngrid,llm),    &  ! diagonal coefficients 
     124                  Ru(ngrid,llm), Rv(ngrid,llm), & ! right-hand-sides 
     125                  C(ngrid,llm),    &  ! LU factorization (Thomas algorithm) 
     126                  xu(ngrid,llm), xv(ngrid, llm)      ! solution 
     127    REAL(rstd) :: rho, X_ij 
     128    INTEGER :: l,ij 
     129    ! Vertical diffusion : rho.du/dt = d/dz (rho*kappa*du/dz)  
     130    ! rho.dz = mass.deta 
     131    ! => mass.du/dt = d/deta ((rho^2 kappa/m)du/deta) 
     132    ! Backward Euler : (M-S)u_new = M.u_old 
     133    ! with M.u = mass(l)*u(l) 
     134    ! and  S.u = A(l+1)*(u(l+1)-u(l)) - A(l)*(u(l)-u(l-1)) 
     135    ! => solve -A(l)u(l-1) + B(l)u(l) - A(l+1)u(l+1) = Ru(l) using Thomas algorithm 
     136    ! with A=tau*kappa*rho^2/m and B(l) = mass(l)+A(l)+A(l+1) 
     137    DO l=1,llm 
     138       !DIR$ SIMD 
     139       DO ij=1,ngrid 
     140          mass(ij,l) = (p(ij,l)-p(ij,l+1))*(1./g) 
     141          rho = (p(ij,l)-p(ij,l+1))/(geopot(ij,l+1)-geopot(ij,l)) 
     142          ! A = kappa.tau.rho^2/m 
     143          A(ij,l) = (kv*dt_phys)*(rho**2)/mass(ij,l) 
     144          Ru(ij,l) = mass(ij,l)*u(ij,l) 
     145          Rv(ij,l) = mass(ij,l)*v(ij,l) 
     146       ENDDO 
     147    ENDDO 
     148    A(:,llm+1)=0. 
     149    DO l=llm,2,-1 
     150       !DIR$ SIMD 
     151       DO ij=1,ngrid 
     152          A(ij,l) = .5*(A(ij,l)+A(ij,l-1)) ! average A to interfaces 
     153       ENDDO 
     154    ENDDO 
     155    A(:,1)=0. 
     156    DO l=1,llm 
     157       !DIR$ SIMD 
     158       DO ij=1,ngrid 
     159          B(ij,l) = mass(ij,l)+A(ij,l)+A(ij,l+1) 
     160       ENDDO 
     161    ENDDO 
     162 
     163    ! Solve -A(l)x(l-1) + B(l)x(l) - A(l+1)x(l+1) = R(l) using Thomas algorithm 
     164    ! Forward sweep : 
     165    ! C(0)=0, C(l) = -A(l+1) / (B(l)+A(l)C(l-1)), 
     166    ! D(0)=0, D(l) = (R(l)+A(l)D(l-1)) / (B(l)+A(l)C(l-1)) 
     167    !DIR$ SIMD 
     168    DO ij=1,ngrid 
     169       X_ij = 1./B(ij,1) 
     170       C(ij,2) = -A(ij,2) * X_ij 
     171       xu(ij,1) = Ru(ij,1) * X_ij 
     172       xv(ij,1) = Rv(ij,1) * X_ij 
     173    END DO 
     174    DO l = 2,llm 
     175       !DIR$ SIMD 
     176       DO ij=1,ngrid 
     177          X_ij = 1./( B(ij,l) + A(ij,l)*C(ij,l) ) 
     178          C(ij,l+1) = -A(ij,l+1) * X_ij ! zero for l=llm 
     179          xu(ij,l) = (Ru(ij,l)+A(ij,l)*xu(ij,l-1)) * X_ij 
     180          xv(ij,l) = (Rv(ij,l)+A(ij,l)*xv(ij,l-1)) * X_ij 
     181       END DO 
     182    END DO 
     183    ! Back substitution : 
     184    ! x(i) = D(i)-C(i+1)x(i+1), x(llm)=D(llm) 
     185    !DIR$ SIMD 
     186    DO ij=1,ngrid 
     187       ! top layer l=llm 
     188       du(ij,llm) = (xu(ij,llm)-u(ij,llm))/dt_phys 
     189       dv(ij,llm) = (xv(ij,llm)-v(ij,llm))/dt_phys 
     190    END DO 
     191    ! Back substitution at lower layers 
     192    DO l = llm-1,1,-1 
     193       !DIR$ SIMD 
     194       DO ij=1,ngrid 
     195          xu(ij,l) = xu(ij,l) - C(ij,l+1)*xu(ij,l+1) 
     196          xv(ij,l) = xv(ij,l) - C(ij,l+1)*xv(ij,l+1) 
     197          du(ij,l) = (xu(ij,l)-u(ij,l))/dt_phys 
     198          dv(ij,l) = (xv(ij,l)-v(ij,l))/dt_phys 
     199       END DO 
     200    END DO 
     201    ! bottom friction + thermal relaxation 
     202    du(:,1)=du(:,1)-kfrict*u(:,1) 
     203    dTemp(:,:) = (1./tauCLee)*(temp_eq_packed(:,:)-temp(:,:)) 
     204  END SUBROUTINE compute_physics_column 
     205 
    87206!----------------------------- Initialize to T_eq -------------------------------------- 
    88207 
    89208  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    90     USE icosa 
    91209    USE theta2theta_rhodz_mod 
    92     IMPLICIT NONE 
    93210    TYPE(t_field),POINTER :: f_ps(:) 
    94211    TYPE(t_field),POINTER :: f_phis(:) 
    … …  
    103220    REAL(rstd),POINTER :: u(:,:) 
    104221    REAL(rstd),POINTER :: q(:,:,:) 
    105     REAL(rstd) :: lat(iim*jjm)        ! latitude                     
    106     REAL(rstd) :: pplay(iim*jjm, llm) ! pressure at full layers 
    107222    INTEGER :: ind 
    108223 
    109     CALL allocate_field(f_temp,field_t,type_real,llm) 
     224    CALL allocate_field(f_temp,field_t,type_real,llm, name='temp_buf_venus') 
    110225    DO ind=1,ndomain 
    111226       IF (.NOT. assigned_domain(ind)) CYCLE 
    … …  
    121236       q(:,:,:)=1e2 
    122237       temp=f_temp(ind) 
    123        CALL compute_temp_ref(temp, .FALSE.) ! Without meridional gradient 
     238       CALL compute_temp_ref(.FALSE., iim*jjm, lat_i, temp) ! Without meridional gradient 
    124239    ENDDO 
    125240    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz) 
    … …  
    130245!------------------------- Compute reference temperature field ------------------------ 
    131246 
    132   SUBROUTINE compute_temp_ref(theta_eq, gradient) 
    133     USE icosa 
     247  SUBROUTINE compute_temp_ref(gradient,ngrid,lat, temp_eq) 
    134248    USE disvert_mod 
    135249    USE omp_para 
    136250    USE math_const 
    137     IMPLICIT NONE 
    138     REAL(rstd), INTENT(OUT) :: theta_eq(iim*jjm,llm)  
    139251    LOGICAL, INTENT(IN) :: gradient 
    140     REAL(rstd)  :: clat(iim*jjm)        ! latitude 
    141     integer :: level 
    142     REAL ::  pressCLee(31), tempCLee(31), dt_epCLee(31), etaCLee(31) 
     252    INTEGER, INTENT(IN) :: ngrid 
     253    REAL(rstd), INTENT(IN) :: lat(ngrid) ! latitude 
     254    REAL(rstd), INTENT(OUT) :: temp_eq(ngrid,llm)  
     255 
     256    REAL(rstd)  :: clat(ngrid)         
     257    INTEGER :: level 
     258    INTEGER, PARAMETER :: nlev=30 
     259    REAL ::  pressCLee(nlev+1), tempCLee(nlev+1), dt_epCLee(nlev+1), etaCLee(nlev+1) 
    143260     
    144     real(rstd) ::  lon,lat, pplay, ztemp,zdt,fact 
    145     logical, save ::  firstcall 
    146     integer :: i,j,ij, l,ll 
     261    REAL(rstd) ::  pplay, ztemp,zdt,fact 
     262    INTEGER :: ij, l,ll 
    147263     
    148264    data etaCLee / 9.602e-1, 8.679e-1, 7.577e-1, 6.420e-1, 5.299e-1, & 
    … …  
    154270                      2.265e-6/ 
    155271 
    156     DO j=jj_begin-1,jj_end+1 
    157        DO i=ii_begin-1,ii_end+1 
    158           ij=(j-1)*iim+i 
    159           CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
    160           clat(ij)=cos(lat)  
    161        ENDDO 
    162     ENDDO 
    163272 
    164273    data   tempCLee/ 728.187, 715.129, 697.876, 677.284, 654.078, 628.885, & 
    … …  
    175284                     2.000/ 
    176285       
    177     pressCLee = etaCLee*9.2e6 
    178      
    179     DO j=jj_begin-1,jj_end+1 
    180        DO i=ii_begin-1,ii_end+1 
    181           ij=(j-1)*iim+i 
    182  
    183           DO l = 1, llm 
    184  
    185              pplay = .5*(ap(l)+ap(l+1)+(bp(l)+bp(l+1))*preff) ! ps=preff           
    186              ! look for largest level such that pressCLee(level) > pplay(ij,l))   
    187              ! => pressClee(level+1) < pplay(ij,l) < pressClee(level) 
    188              level = 1 
    189              DO ll = 1, 30 ! 30 data levels 
    190                 IF(pressCLee(ll) > pplay) THEN 
    191                    level = ll 
    192                 END IF 
    193              END DO 
    194               
    195              ! interpolate between level and level+1 
    196              ! interpolation is linear in log(pressure) 
    197              fact  = ( log10(pplay)-log10(pressCLee(level))) & 
    198                     /( log10(pressCLee(level+1))-log10(pressCLee(level)) ) 
    199              ztemp = tempCLee(level)*(1-fact) + tempCLee(level+1)*fact 
    200              zdt   = dt_epCLee(level)*(1-fact) + dt_epCLee(level+1)*fact 
    201  
    202              IF(gradient) THEN 
    203                 theta_eq(ij,l) = ztemp+ zdt*(clat(ij)-Pi/4.) 
    204              ELSE 
    205                 theta_eq(ij,l) = ztemp 
     286    pressCLee(:) = etaCLee(:)*9.2e6 
     287    clat(:)=COS(lat(:)) 
     288 
     289    DO ij=1,ngrid 
     290       DO l = 1, llm 
     291         
     292          pplay = .5*(ap(l)+ap(l+1)+(bp(l)+bp(l+1))*preff) ! ps=preff           
     293          ! look for largest level such that pressCLee(level) > pplay(ij,l))   
     294          ! => pressClee(level+1) < pplay(ij,l) < pressClee(level) 
     295          level = 1 
     296          DO ll = 1, nlev ! nlev data levels 
     297             IF(pressCLee(ll) > pplay) THEN 
     298                level = ll 
    206299             END IF 
    207300          END DO 
     301           
     302          ! interpolate between level and level+1 
     303          ! interpolation is linear in log(pressure) 
     304          fact  = ( log10(pplay)-log10(pressCLee(level))) & 
     305               /( log10(pressCLee(level+1))-log10(pressCLee(level)) ) 
     306          ztemp = tempCLee(level)*(1-fact) + tempCLee(level+1)*fact 
     307          zdt   = dt_epCLee(level)*(1-fact) + dt_epCLee(level+1)*fact 
     308           
     309          IF(gradient) THEN 
     310             temp_eq(ij,l) = ztemp+ zdt*(clat(ij)-Pi/4.) 
     311          ELSE 
     312             temp_eq(ij,l) = ztemp 
     313          END IF 
    208314       END DO 
    209315    END DO 

  • codes/icosagcm/devel/src/physics/physics.f90

    r739 r741  
    1212       phys_lmdz_generic=21, phys_external=22 
    1313  INTEGER :: phys_type 
    14   TYPE(t_field),POINTER,SAVE :: f_extra_physics_2D(:), f_extra_physics_3D(:) 
    1514  TYPE(t_field),POINTER,SAVE :: f_dulon(:), f_dulat(:) 
    1615  TYPE(t_field),POINTER,SAVE :: f_ulon(:), f_ulat(:) 
    … …  
    5049    CASE ('held_suarez') 
    5150       phys_type = phys_HS94 
    52     CASE ('Lebonnois2012') 
    53        phys_type = phys_LB2012 
    54        CALL init_phys_venus        
    5551    CASE ('phys_lmdz_generic') 
    5652       phys_type=phys_lmdz_generic 
    … …  
    7975          phys_type = phys_DCMIP2016 
    8076          CALL init_physics_dcmip2016 
     77       CASE ('Lebonnois2012') 
     78          phys_type = phys_LB2012 
     79          CALL init_phys_venus        
    8180       CASE DEFAULT 
    8281          IF(is_mpi_root) PRINT*, 'init_physics : Bad selector for variable physics <',& 
    … …  
    129128  END SUBROUTINE add_du_phys 
    130129 
    131   SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q) 
     130  SUBROUTINE physics(it,f_phis, f_geopot, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q) 
    132131    USE physics_lmdz_generic_mod, ONLY : physics_lmdz_generic => physics 
    133132    USE physics_external_mod, ONLY : physics_external => physics 
    … …  
    135134    USE physics_dcmip2016_mod, ONLY : write_physics_dcmip2016 => write_physics 
    136135    USE etat0_heldsz_mod 
    137     USE etat0_venus_mod, ONLY : phys_venus => physics 
    138136    INTEGER, INTENT(IN)   :: it 
    139137    TYPE(t_field),POINTER :: f_phis(:) 
     138    TYPE(t_field),POINTER :: f_geopot(:) 
    140139    TYPE(t_field),POINTER :: f_ps(:) 
    141140    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    … …  
    161160       CASE (phys_external) 
    162161         CALL physics_external(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q) 
    163        CASE(phys_LB2012) 
    164           CALL phys_venus(f_ps,f_theta_rhodz,f_ue)  
    165162       CASE DEFAULT 
    166           CALL physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q) 
     163          CALL physics_column(it, f_phis, f_geopot, f_ps, f_theta_rhodz, f_ue, f_q) 
    167164       END SELECT 
    168165 
    … …  
    197194  END SUBROUTINE write_physics_tendencies 
    198195     
    199   SUBROUTINE physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q) 
     196  SUBROUTINE physics_column(it, f_phis, f_geopot, f_ps, f_theta_rhodz, f_ue, f_q) 
    200197    USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics 
    201198    USE physics_dcmip2016_mod, ONLY : full_physics_dcmip2016 => full_physics 
     199    USE etat0_venus_mod, ONLY : full_physics_venus=>full_physics 
    202200    USE theta2theta_rhodz_mod 
    203201    USE mpipara 
    204202    USE checksum_mod 
    205203    TYPE(t_field),POINTER :: f_phis(:) 
     204    TYPE(t_field),POINTER :: f_geopot(:) 
    206205    TYPE(t_field),POINTER :: f_ps(:) 
    207206    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    … …  
    209208    TYPE(t_field),POINTER :: f_q(:) 
    210209    REAL(rstd),POINTER :: phis(:) 
     210    REAL(rstd),POINTER :: geopot(:,:) 
    211211    REAL(rstd),POINTER :: ps(:) 
    212212    REAL(rstd),POINTER :: temp(:,:) 
    … …  
    228228       CALL swap_geometry(ind) 
    229229       phis=f_phis(ind) 
     230       geopot=f_geopot(ind) 
    230231       ps=f_ps(ind) 
    231232       temp=f_temp(ind) 
    … …  
    236237       ulon=f_ulon(ind) 
    237238       ulat=f_ulat(ind) 
    238        CALL pack_physics(pack_info(ind), phis, ps, temp, ue, q, p, pk, ulon, ulat) 
     239       CALL pack_physics(pack_info(ind), phis, geopot, ps, temp, ue, q, p, pk, ulon, ulat) 
    239240    END DO 
    240241 
    … …  
    244245    CASE (phys_DCMIP2016) 
    245246       IF (is_omp_level_master) CALL full_physics_dcmip2016 
     247    CASE(phys_LB2012) 
     248       IF (is_omp_level_master) CALL full_physics_venus 
    246249    CASE DEFAULT 
    247250       IF(is_master) PRINT *,'Internal error : illegal value of phys_type', phys_type 
    … …  
    279282  END SUBROUTINE physics_column 
    280283 
    281   SUBROUTINE pack_physics(info, phis, ps, temp, ue, q, p, pk, ulon, ulat ) 
     284  SUBROUTINE pack_physics(info, phis, geopot, ps, temp, ue, q, p, pk, ulon, ulat ) 
    282285    USE wind_mod 
    283286    USE pression_mod 
    … …  
    286289    TYPE(t_pack_info) :: info 
    287290    REAL(rstd) :: phis(iim*jjm) 
     291    REAL(rstd) :: geopot(iim*jjm,llm+1) 
    288292    REAL(rstd) :: ps(iim*jjm) 
    289293    REAL(rstd) :: temp(iim*jjm,llm) 
    … …  
    308312    IF (is_omp_level_master) THEN 
    309313      CALL pack_domain(info, phis, physics_inout%phis) 
     314      CALL pack_domain(info, geopot, physics_inout%geopot) 
    310315      CALL pack_domain(info, p, physics_inout%p) 
    311316      CALL pack_domain(info, pk, physics_inout%pk) 

  • codes/icosagcm/devel/src/physics/physics_interface.f90

    r739 r741  
    1111     REAL(rstd), DIMENSION(:), POINTER :: Ai, lon, lat, phis 
    1212     ! Input, time-dependent 
    13      REAL(rstd), DIMENSION(:,:), POINTER :: p, pk, Temp, ulon, ulat 
     13     REAL(rstd), DIMENSION(:,:), POINTER :: geopot, p, pk, Temp, ulon, ulat 
    1414     REAL(rstd), DIMENSION(:,:,:), POINTER :: q 
    1515     ! Output arrays 
    … …  
    9393 
    9494    ngrid=offset 
    95     ! Input                                                                                                      
     95    ! Input 
    9696    ALLOCATE(physics_inout%Ai(ngrid)) 
    9797    ALLOCATE(physics_inout%lon(ngrid)) 
    … …  
    9999    ALLOCATE(physics_inout%phis(ngrid)) 
    100100    ALLOCATE(physics_inout%p(ngrid,llm+1)) 
     101    ALLOCATE(physics_inout%geopot(ngrid,llm+1)) 
    101102    ALLOCATE(physics_inout%pk(ngrid,llm)) 
    102103    ALLOCATE(physics_inout%Temp(ngrid,llm)) 

  • codes/icosagcm/devel/src/time/timeloop_gcm.f90

    r732 r741  
    350350            f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 
    351351          CALL enter_profile(id_phys) 
    352           CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)         
     352          CALL physics(it, f_phis, f_geopot, f_ps, f_theta_rhodz, f_u, f_wflux, f_q) 
    353353          CALL exit_profile(id_phys) 
    354354          CALL check_conserve_detailed(it, AAM_phys, & 
Note: See TracChangeset for help on using the changeset viewer.