Changeset 167


Ignore:
Timestamp:
06/29/13 02:35:03 (11 years ago)
Author:
dubos
Message:

Multi-layer Saint-Venant / Ripa equations - tested with Williamson91.5 & JW06

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

Legend:

Unmodified
Added
Removed

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

    r162 r167  
    2222 
    2323  PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields, & 
    24        req_ps, req_mass, caldyn_eta, eta_mass, eta_lag 
     24       req_ps, req_mass 
    2525 
    2626CONTAINS 
    … …  
    122122       f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du) 
    123123    USE icosa 
    124     USE disvert_mod 
     124    USE disvert_mod, ONLY : caldyn_eta, eta_mass 
    125125    USE vorticity_mod 
    126126    USE kinetic_mod 
    … …  
    285285SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv) 
    286286  USE icosa 
    287   USE disvert_mod 
     287  USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass 
    288288  USE exner_mod 
    289289  USE trace 
    … …  
    318318           DO i=ii_begin-1,ii_end+1 
    319319              ij=(j-1)*iim+i 
    320               m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g  
    321 !              m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 
     320              m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 
    322321              rhodz(ij,l) = m 
    323322              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    … …  
    420419       ENDDO 
    421420 
    422     ELSE 
    423  
     421    ELSE  
    424422       ! We are using a Lagrangian vertical coordinate 
    425423       ! Pressure must be computed first top-down (temporarily stored in pk) 
    426424       ! Then Exner pressure and geopotential are computed bottom-up 
    427        ! Notice that the computation below should work also when caldyn_eta=eta_mass 
     425       ! Notice that the computation below should work also when caldyn_eta=eta_mass           
    428426 
    429427       ! uppermost layer 
    … …  
    451449          END DO 
    452450       END DO 
    453        DO l = 1,llm 
    454           !$OMP DO SCHEDULE(STATIC)  
    455           DO j=jj_begin-1,jj_end+1 
    456              DO i=ii_begin-1,ii_end+1 
    457                 ij=(j-1)*iim+i 
    458                 p_ik = pk(ij,l) 
    459                 exner_ik = cpp * (p_ik/preff) ** kappa 
    460                 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
    461                 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik  
    462                 pk(ij,l) = exner_ik 
     451 
     452       IF(boussinesq) THEN ! specific volume 1 = dphi/g/rhodz 
     453          DO l = 1,llm 
     454             !$OMP DO SCHEDULE(STATIC)  
     455             DO j=jj_begin-1,jj_end+1 
     456                DO i=ii_begin-1,ii_end+1 
     457                   ij=(j-1)*iim+i 
     458                   geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l) 
     459                ENDDO 
    463460             ENDDO 
    464461          ENDDO 
    465        ENDDO 
     462       ELSE ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
     463          DO l = 1,llm 
     464             !$OMP DO SCHEDULE(STATIC)  
     465             DO j=jj_begin-1,jj_end+1 
     466                DO i=ii_begin-1,ii_end+1 
     467                   ij=(j-1)*iim+i 
     468                   p_ik = pk(ij,l) 
     469                   exner_ik = cpp * (p_ik/preff) ** kappa 
     470                   geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik  
     471                   pk(ij,l) = exner_ik 
     472                ENDDO 
     473             ENDDO 
     474          ENDDO 
     475       END IF 
    466476 
    467477    END IF 
    … …  
    482492    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm) 
    483493    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature 
    484     REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm) ! Exner function 
     494    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function 
    485495    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential 
    486496 
    … …  
    647657   
    648658!!! Compute bernouilli term = Kinetic Energy + geopotential 
    649    DO l=ll_begin,ll_end 
    650     DO j=jj_begin,jj_end 
    651       DO i=ii_begin,ii_end 
    652         ij=(j-1)*iim+i 
    653  
    654         berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
     659    IF(boussinesq) THEN 
     660 
     661       DO l=ll_begin,ll_end 
     662          DO j=jj_begin,jj_end 
     663             DO i=ii_begin,ii_end 
     664                ij=(j-1)*iim+i 
     665                 
     666                berni(ij,l) = pk(ij,l) + & 
     667                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    & 
     668                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
     669                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
     670                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       & 
     671                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
     672                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
     673                pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 
     674             ENDDO 
     675          ENDDO 
     676       ENDDO 
     677 
     678    ELSE 
     679 
     680       DO l=ll_begin,ll_end 
     681          DO j=jj_begin,jj_end 
     682             DO i=ii_begin,ii_end 
     683                ij=(j-1)*iim+i 
     684                 
     685                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
    655686                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    & 
    656687                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
    … …  
    660691                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    661692        
    662       ENDDO 
    663     ENDDO 
    664   ENDDO 
    665  
    666   
     693             ENDDO 
     694          ENDDO 
     695       ENDDO 
     696 
     697    END IF ! Boussinesq/compressible 
     698 
    667699!!! Add gradients of Bernoulli and Exner functions to du 
    668  
    669    DO l=ll_begin,ll_end 
    670     DO j=jj_begin,jj_end 
    671       DO i=ii_begin,ii_end 
    672         ij=(j-1)*iim+i 
     700    DO l=ll_begin,ll_end 
     701       DO j=jj_begin,jj_end 
     702          DO i=ii_begin,ii_end 
     703             ij=(j-1)*iim+i 
    673704         
    674         du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             & 
    675                           0.5*(theta(ij,l)+theta(ij+t_right,l))                & 
    676                              *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    & 
    677                         + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) ) 
     705             du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             & 
     706                               0.5*(theta(ij,l)+theta(ij+t_right,l))                & 
     707                                  *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    & 
     708                                   + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) ) 
    678709         
    679710         
    680         du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  & 
    681                           0.5*(theta(ij,l)+theta(ij+t_lup,l))                  & 
    682                              *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       & 
    683                         + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) ) 
     711             du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  & 
     712                               0.5*(theta(ij,l)+theta(ij+t_lup,l))                  & 
     713                                  *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       & 
     714                                   + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) ) 
    684715                 
    685         du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             & 
    686                           0.5*(theta(ij,l)+theta(ij+t_ldown,l))                & 
    687                              *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     & 
    688                         + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) ) 
    689  
    690       ENDDO 
     716             du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             & 
     717                               0.5*(theta(ij,l)+theta(ij+t_ldown,l))                & 
     718                                  *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     & 
     719                                   + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) ) 
     720 
     721          ENDDO 
     722       ENDDO 
    691723    ENDDO 
    692   ENDDO 
    693  
    694   CALL trace_end("compute_caldyn_horiz") 
     724  
     725    CALL trace_end("compute_caldyn_horiz") 
    695726 
    696727END SUBROUTINE compute_caldyn_horiz 

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

    r151 r167  
    2727  REAL, save    :: rayleigh_tau 
    2828   
    29 !  INTEGER,SAVE :: itau_dissip 
    3029  REAL,SAVE    :: dtdissip 
    3130  
    … …  
    368367    fact=2 
    369368    DO l=1,llm 
    370        zz= 1. - preff/presnivs(l) 
     369       IF(ap_bp_present) THEN ! height-dependent dissipation 
     370          zz= 1. - preff/presnivs(l) 
     371       ELSE 
     372          zz = 0. 
     373       END IF 
    371374       zvert=fact-(fact-1)/(1+zz*zz) 
    372375       tau_graddiv(l) = tau_graddiv(l)/zvert 
    … …  
    390393  
    391394   
    392   SUBROUTINE dissip(f_ue,f_due,f_ps,f_phis,f_theta_rhodz,f_dtheta_rhodz) 
     395  SUBROUTINE dissip(f_ue,f_due,f_mass,f_phis,f_theta_rhodz,f_dtheta_rhodz) 
    393396  USE icosa 
    394397  USE theta2theta_rhodz_mod 
    … …  
    402405    TYPE(t_field),POINTER :: f_ue(:) 
    403406    TYPE(t_field),POINTER :: f_due(:) 
    404     TYPE(t_field),POINTER :: f_ps(:), f_phis(:) 
     407    TYPE(t_field),POINTER :: f_mass(:), f_phis(:) 
    405408    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    406409    TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 
    … …  
    421424    CALL gradiv(f_ue,f_due_diss1) 
    422425    CALL gradrot(f_ue,f_due_diss2) 
    423      
    424     CALL divgrad_theta_rhodz(f_ps,f_theta_rhodz,f_dtheta_rhodz_diss) 
     426 
     427    CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss) 
    425428     
    426429! later for openmp     
    … …  
    646649  END SUBROUTINE divgrad 
    647650    
    648   SUBROUTINE divgrad_theta_rhodz(f_ps,f_theta_rhodz,f_dtheta_rhodz) 
     651  SUBROUTINE divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz) 
    649652  USE icosa 
    650653  USE trace 
    651654  USE omp_para 
    652   USE disvert_mod 
    653655  IMPLICIT NONE 
    654     TYPE(t_field),POINTER :: f_ps(:) 
     656    TYPE(t_field),POINTER :: f_mass(:) 
    655657    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    656658    TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 
    657659     
    658     REAL(rstd),POINTER :: ps(:) 
     660    REAL(rstd),POINTER :: mass(:,:) 
    659661    REAL(rstd),POINTER :: theta_rhodz(:,:) 
    660662    REAL(rstd),POINTER :: dtheta_rhodz(:,:) 
    … …  
    668670      CALL swap_dimensions(ind) 
    669671      CALL swap_geometry(ind) 
    670       ps=f_ps(ind) 
     672      mass=f_mass(ind) 
    671673      theta_rhodz=f_theta_rhodz(ind) 
    672674      dtheta_rhodz=f_dtheta_rhodz(ind)  
    … …  
    675677          DO i=ii_begin,ii_end 
    676678            ij=(j-1)*iim+i 
    677             dtheta_rhodz(ij,l) = theta_rhodz(ij,l) / ( ap(l)-ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g 
     679            dtheta_rhodz(ij,l) = theta_rhodz(ij,l) / mass(ij,l) 
    678680          ENDDO 
    679681        ENDDO 
    … …  
    703705          DO i=ii_begin,ii_end 
    704706            ij=(j-1)*iim+i 
    705             dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * ( ap(l)-ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g 
     707            dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * mass(ij,l) 
    706708          ENDDO 
    707709        ENDDO 

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

    r162 r167  
    9696!---------------------------------------------------- 
    9797 
    98  
    99     def='eta_mass' 
    100     CALL getin('caldyn_eta',def) 
    101     SELECT CASE(TRIM(def)) 
    102     CASE('eta_mass') 
    103        caldyn_eta=eta_mass 
    104     CASE('eta_lag') 
    105        caldyn_eta=eta_lag 
    106     CASE DEFAULT 
    107        IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_eta : <', TRIM(def),'> options are <eta_mass>, <eta_lag>' 
    108        STOP 
    109     END SELECT 
    110     IF (is_mpi_root) PRINT *, 'caldyn_eta=',TRIM(def), caldyn_eta, eta_lag, eta_mass 
    11198 
    11299! Time-independant orography 
    … …  
    284271 
    285272    IF (MOD(it+1,itau_dissip)==0) THEN 
    286       CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz) 
    287       CALL euler_scheme(.FALSE.) 
    288     ENDIF 
     273       IF(caldyn_eta==eta_mass) THEN 
     274          DO ind=1,ndomain 
     275             CALL swap_dimensions(ind) 
     276             CALL swap_geometry(ind) 
     277             mass=f_mass(ind); ps=f_ps(ind); 
     278             CALL compute_rhodz(.TRUE., ps, mass) 
     279          END DO 
     280       ENDIF 
     281       CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz) 
     282       CALL euler_scheme(.FALSE.)  ! update only u, theta 
     283    END IF 
    289284 
    290285    IF(MOD(it+1,itau_adv)==0) THEN 
Note: See TracChangeset for help on using the changeset viewer.