Changeset 156


Ignore:
Timestamp:
06/03/13 23:45:04 (11 years ago)
Author:
dubos
Message:

new hydrostatic balance - tested with JW06 - code cleanup to follow

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

Legend:

Unmodified
Added
Removed

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

    r151 r156  
    1414  TYPE(t_field),POINTER :: f_theta(:) 
    1515  TYPE(t_field),POINTER :: f_pk(:) 
    16   TYPE(t_field),POINTER :: f_pks(:) 
     16!  TYPE(t_field),POINTER :: f_pks(:) 
    1717  TYPE(t_field),POINTER :: f_phi(:) 
     18  TYPE(t_field),POINTER :: f_geopot(:) 
    1819  TYPE(t_field),POINTER :: f_divm(:) 
    1920  TYPE(t_field),POINTER :: f_wwuu(:) 
    … …  
    9495    CALL allocate_field(f_theta,field_t,type_real,llm)  ! potential temperature 
    9596    CALL allocate_field(f_pk,field_t,type_real,llm) 
    96     CALL allocate_field(f_pks,field_t,type_real) ! Exner function 
    97     CALL allocate_field(f_phi,field_t,type_real,llm)    ! geopotential 
     97!    CALL allocate_field(f_pks,field_t,type_real) ! Exner function 
     98!    CALL allocate_field(f_phi,field_t,type_real,llm)     ! geopotential 
     99    CALL allocate_field(f_geopot,field_t,type_real,llm+1)  ! geopotential (new) 
    98100    CALL allocate_field(f_divm,field_t,type_real,llm)  ! mass flux divergence 
    99101    CALL allocate_field(f_wwuu,field_u,type_real,llm+1)  
    … …  
    129131! temporary shared variable 
    130132    REAL(rstd),POINTER  :: theta(:,:)   
    131     REAL(rstd),POINTER  :: pk(:,:), pks(:) 
     133    REAL(rstd),POINTER  :: pk(:,:) !, pks(:) 
    132134    REAL(rstd),POINTER  :: phi(:,:)     
     135    REAL(rstd),POINTER  :: geopot(:,:) 
    133136    REAL(rstd),POINTER  :: divm(:,:)  
    134137    REAL(rstd),POINTER  :: wwuu(:,:) 
    … …  
    189192          theta = f_theta(ind) 
    190193          pk = f_pk(ind) 
    191           pks = f_pks(ind) 
    192           phi = f_phi(ind)   
     194!          pks = f_pks(ind) 
     195!          phi = f_phi(ind) 
     196          geopot = f_geopot(ind)   
    193197          divm = f_divm(ind) 
    194198          wwuu=f_wwuu(ind) 
    195199 
    196200          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du, & 
    197                               theta,pk, pks, phi, divm, wwuu) 
     201                              theta,pk, geopot, divm, wwuu) 
     202!                              theta,pk, pks, phi, geopot, divm, wwuu) 
    198203       ENDDO        
    199204        
    … …  
    203208          CALL swap_geometry(ind) 
    204209          phis=f_phis(ind) 
     210          phi=f_phi(ind) 
     211          geopot = f_geopot(ind)   
    205212          ps=f_ps(ind) 
    206213          dps=f_dps(ind) 
    … …  
    218225          CALL compute_pvort(ps, u, p,rhodz,qu) 
    219226          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du,   & 
    220                               theta,pk, pks, phi, divm, wwuu) 
     227                              theta,pk, phi, divm, wwuu) 
     228!                              theta,pk, pks, phi, geopot, divm, wwuu) 
    221229       ENDDO 
    222230        
    … …  
    342350   
    343351  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du,   & 
    344                             theta,pk, pks, phi, divm, wwuu)  
     352                            theta,pk, geopot, divm, wwuu)  
     353!                            theta,pk, pks, phi, geopot, divm, wwuu)  
    345354  USE icosa 
    346355  USE disvert_mod 
    … …  
    364373! temporary variable 
    365374    REAL(rstd),INTENT(INOUT) :: theta(iim*jjm,llm)  ! potential temperature 
    366     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm), pks(iim*jjm) ! Exner function 
    367     REAL(rstd),INTENT(INOUT) :: phi(iim*jjm,llm)    ! geopotential 
     375    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) !, pks(iim*jjm) ! Exner function 
     376!    REAL(rstd),INTENT(INOUT) :: phi(iim*jjm,llm)    ! geopotential 
     377    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)    ! geopotential 
    368378    REAL(rstd),INTENT(INOUT) :: divm(iim*jjm,llm)  ! mass flux divergence 
    369379    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1) 
    370  
    371380     
    372381    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux 
    … …  
    374383     
    375384    INTEGER :: i,j,ij,l 
    376     REAL(rstd) :: ww,uu 
     385    REAL(rstd) :: ww,uu, p_ik, exner_ik 
    377386 
    378387    CALL trace_start("compute_caldyn") 
    … …  
    587596    END SELECT 
    588597   
    589 !!! Compute Exner function 
    590 !    CALL compute_exner(ps,p,pks,pk,1)  
    591 ! replaced in source 
    592     IF (omp_first) THEN     
     598!!!! Compute Exner function 
     599!!    CALL compute_exner(ps,p,pks,pk,1)  
     600!! replaced in source 
     601!    IF (omp_first) THEN     
     602!       DO j=jj_begin-1,jj_end+1 
     603!          DO i=ii_begin-1,ii_end+1 
     604!             ij=(j-1)*iim+i 
     605!             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 
     606!          ENDDO 
     607!       ENDDO 
     608!     ENDIF 
     609!      
     610!       ! 3D : pk 
     611!     DO l = ll_begin, ll_end 
     612!        DO j=jj_begin-1,jj_end+1 
     613!           DO i=ii_begin-1,ii_end+1 
     614!              ij=(j-1)*iim+i 
     615!              pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa 
     616!           ENDDO 
     617!        ENDDO 
     618!     ENDDO 
     619! 
     620!!---> flush pk,pks, theta 
     621!!$OMP BARRIER 
     622! 
     623!!! Compute geopotential (old) 
     624! 
     625!  ! for first layer 
     626!  IF (omp_first) THEN 
     627!    DO j=jj_begin-1,jj_end+1 
     628!      DO i=ii_begin-1,ii_end+1 
     629!        ij=(j-1)*iim+i 
     630!        phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) ) 
     631!      ENDDO 
     632!    ENDDO 
     633!  ENDIF 
     634!!!-> implicit flush on phi(:,1) 
     635!           
     636!  ! for other layers 
     637!   DO l = ll_beginp1, ll_end 
     638!     DO j=jj_begin-1,jj_end+1 
     639!       DO i=ii_begin-1,ii_end+1 
     640!         ij=(j-1)*iim+i 
     641!         phi(ij,l) =  0.5 * ( theta(ij,l)  + theta(ij,l-1) )  &  
     642!                          * (  pk(ij,l-1) -  pk(ij,l)    ) 
     643!       ENDDO 
     644!     ENDDO 
     645!   ENDDO 
     646! 
     647!!$OMP BARRIER 
     648!   DO l = 2, llm 
     649!!$OMP DO 
     650!     DO j=jj_begin-1,jj_end+1 
     651!       DO i=ii_begin-1,ii_end+1 
     652!         ij=(j-1)*iim+i 
     653!         phi(ij,l) = phi(ij,l)+ phi(ij,l-1) 
     654!       ENDDO 
     655!     ENDDO 
     656!   ENDDO 
     657!! --> IMPLICIT FLUSH on phi 
     658       
     659!!! Compute exner function and geopotential 
     660 
     661! geopot=phis for first layer 
     662!$OMP DO SCHEDULE(STATIC)  
     663    DO j=jj_begin-1,jj_end+1 
     664       DO i=ii_begin-1,ii_end+1 
     665          ij=(j-1)*iim+i 
     666          geopot(ij,1) = phis(ij) 
     667       ENDDO 
     668    ENDDO 
     669! for other layers 
     670    DO l = 1,llm 
     671!$OMP DO SCHEDULE(STATIC)  
    593672       DO j=jj_begin-1,jj_end+1 
    594673          DO i=ii_begin-1,ii_end+1 
    595674             ij=(j-1)*iim+i 
    596              pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 
     675             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 
     676             !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j)) 
     677             exner_ik = cpp * (p_ik/preff) ** kappa 
     678             pk(ij,l) = exner_ik 
     679             ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
     680             geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik  
    597681          ENDDO 
    598682       ENDDO 
    599      ENDIF 
    600       
    601        ! 3D : pk 
    602      DO l = ll_begin, ll_end 
    603         DO j=jj_begin-1,jj_end+1 
    604            DO i=ii_begin-1,ii_end+1 
    605               ij=(j-1)*iim+i 
    606               pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa 
    607            ENDDO 
    608         ENDDO 
    609      ENDDO 
    610  
    611 !---> flush pk,pks, theta 
    612 !$OMP BARRIER 
    613  
    614 !!! Compute geopotential 
    615  
    616   ! for first layer 
    617   IF (omp_first) THEN 
    618     DO j=jj_begin-1,jj_end+1 
    619       DO i=ii_begin-1,ii_end+1 
    620         ij=(j-1)*iim+i 
    621         phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) ) 
    622       ENDDO 
    623     ENDDO 
    624   ENDIF 
    625 !!-> implicit flush on phi(:,1) 
    626            
    627   ! for other layers 
    628    DO l = ll_beginp1, ll_end 
    629      DO j=jj_begin-1,jj_end+1 
    630        DO i=ii_begin-1,ii_end+1 
    631          ij=(j-1)*iim+i 
    632          phi(ij,l) =  0.5 * ( theta(ij,l)  + theta(ij,l-1) )  &  
    633                           * (  pk(ij,l-1) -  pk(ij,l)    ) 
    634        ENDDO 
    635      ENDDO 
    636    ENDDO 
    637  
    638 !$OMP BARRIER 
    639    DO l = 2, llm 
    640 !$OMP DO 
    641      DO j=jj_begin-1,jj_end+1 
    642        DO i=ii_begin-1,ii_end+1 
    643          ij=(j-1)*iim+i 
    644          phi(ij,l) = phi(ij,l)+ phi(ij,l-1) 
    645        ENDDO 
    646      ENDDO 
    647    ENDDO 
    648 ! --> IMPLICIT FLUSH on phi 
    649        
     683    ENDDO 
     684 
    650685!!! Compute bernouilli term = Kinetic Energy + geopotential 
    651686   DO l=ll_begin,ll_end 
    … …  
    654689        ij=(j-1)*iim+i 
    655690 
    656         berni(ij,l) =  phi(ij,l) &                                                           
     691        berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
    657692                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    & 
    658693                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 

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

    r131 r156  
    44  REAL(rstd), SAVE, POINTER :: bp(:) 
    55  REAL(rstd), SAVE, POINTER :: presnivs(:) 
     6  REAL(rstd), SAVE, POINTER :: mass_al(:), mass_bl(:), mass_ak(:), mass_bk(:), mass_dak(:), mass_dbk(:) 
     7  REAL(rstd) :: ptop ! pressure at top of atmosphere l=llm+1 
    68 
    79CONTAINS 
    … …  
    1719  IMPLICIT NONE 
    1820    CHARACTER(LEN=255) :: disvert_type = 'std' 
    19      
     21    INTEGER :: l 
     22 
    2023    CALL getin("disvert",disvert_type) 
    2124     
    … …  
    6164         
    6265    END SELECT 
     66 
     67    ! Convert from pressure coordinate to mass coordinate 
     68    ! pk = ap(k) + bp(k)*ps(ij) = ptop + g*( mass_ak(k) + mass_bk(k)*ms(ij) ) 
     69    ptop = ap(llm+1) 
     70    ALLOCATE(mass_al(llm+1)) 
     71    ALLOCATE(mass_bl(llm+1)) 
     72    ALLOCATE(mass_ak(llm)) 
     73    ALLOCATE(mass_bk(llm)) 
     74    ALLOCATE(mass_dak(llm)) 
     75    ALLOCATE(mass_dbk(llm)) 
     76 
     77    ! FIXME : leave ps for the moment ; change ps to ms later 
     78    DO l=1,llm+1 
     79!       mass_al(l) = (ap(l)-ptop)/g 
     80!       mass_bl(l) = bp(l)/g 
     81       mass_al(l) = (ap(l)-ptop) 
     82       mass_bl(l) = bp(l) 
     83    END DO 
     84    DO l=1,llm 
     85       mass_ak(l) = .5*(mass_al(l)+mass_al(l+1)) 
     86       mass_bk(l) = .5*(mass_bl(l)+mass_bl(l+1)) 
     87       mass_dak(l) = mass_al(l)-mass_al(l+1)  ! >0 
     88       mass_dbk(l) = mass_bl(l)-mass_bl(l+1)  ! >0 
     89    END DO 
    6390 
    6491  END SUBROUTINE init_disvert   
Note: See TracChangeset for help on using the changeset viewer.