Changeset 126


Ignore:
Timestamp:
02/01/13 02:11:24 (11 years ago)
Author:
dubos
Message:

Ready for energy-conserving variant

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r125 r126  
    2020    CHARACTER(len=255) :: def 
    2121   
    22     def='direct' 
     22    def='enstrophy' 
    2323    CALL getin('caldyn_conserv',def) 
    2424    SELECT CASE(TRIM(def)) 
    2525    CASE('energy') 
    2626       caldyn_conserv=1 
    27     CASE('direct') 
     27    CASE('enstrophy') 
    2828       caldyn_conserv=2 
    2929    CASE DEFAULT 
    … …  
    7979  END SUBROUTINE allocate_caldyn 
    8080    
    81   SUBROUTINE check_mass_conservation(f_ps,f_dps) 
    82   USE icosa 
    83   IMPLICIT NONE 
     81  SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
     82    USE icosa 
     83    USE vorticity_mod 
     84    USE kinetic_mod 
     85    USE theta2theta_rhodz_mod 
     86    IMPLICIT NONE 
     87    INTEGER,INTENT(IN)    :: it 
     88    TYPE(t_field),POINTER :: f_phis(:) 
    8489    TYPE(t_field),POINTER :: f_ps(:) 
     90    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
     91    TYPE(t_field),POINTER :: f_u(:) 
     92    TYPE(t_field),POINTER :: f_q(:) 
    8593    TYPE(t_field),POINTER :: f_dps(:) 
    86     REAL(rstd),POINTER :: ps(:) 
    87     REAL(rstd),POINTER :: dps(:) 
    88     REAL(rstd) :: mass_tot,dmass_tot 
    89     INTEGER :: ind,i,j,ij 
    90      
    91     mass_tot=0 
    92     dmass_tot=0 
    93      
    94     CALL transfert_request(f_dps,req_i1) 
    95     CALL transfert_request(f_ps,req_i1) 
    96  
    97     DO ind=1,ndomain 
    98       CALL swap_dimensions(ind) 
    99       CALL swap_geometry(ind) 
    100  
    101       ps=f_ps(ind) 
    102       dps=f_dps(ind) 
    103  
    104       DO j=jj_begin,jj_end 
    105         DO i=ii_begin,ii_end 
    106           ij=(j-1)*iim+i 
    107           IF (domain(ind)%own(i,j)) THEN 
    108             mass_tot=mass_tot+ps(ij)*Ai(ij)/g 
    109             dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g 
    110           ENDIF 
    111         ENDDO 
    112       ENDDO 
    113     
    114     ENDDO 
    115     PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot         
    116  
    117   END SUBROUTINE check_mass_conservation 
    118    
    119   SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
    120   USE icosa 
    121   USE vorticity_mod 
    122   USE kinetic_mod 
    123   USE theta2theta_rhodz_mod 
    124   IMPLICIT NONE 
    125   INTEGER,INTENT(IN)    :: it 
    126   TYPE(t_field),POINTER :: f_phis(:) 
    127   TYPE(t_field),POINTER :: f_ps(:) 
    128   TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    129   TYPE(t_field),POINTER :: f_u(:) 
    130   TYPE(t_field),POINTER :: f_q(:) 
    131   TYPE(t_field),POINTER :: f_dps(:) 
    132   TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 
    133   TYPE(t_field),POINTER :: f_du(:) 
    134  
    135   REAL(rstd),POINTER :: phis(:), ps(:), dps(:) 
    136   REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:) 
    137   REAL(rstd),POINTER :: u(:,:), du(:,:) 
    138   REAL(rstd),POINTER :: p(:,:), rhodz(:,:), qu(:,:) 
    139   INTEGER :: ind,ij 
    140  
    141    
     94    TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 
     95    TYPE(t_field),POINTER :: f_du(:) 
     96     
     97    REAL(rstd),POINTER :: phis(:), ps(:), dps(:) 
     98    REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:) 
     99    REAL(rstd),POINTER :: u(:,:), du(:,:) 
     100    REAL(rstd),POINTER :: p(:,:), rhodz(:,:), qu(:,:) 
     101    INTEGER :: ind,ij 
     102     
    142103    CALL transfert_request(f_phis,req_i1)  
    143104    CALL transfert_request(f_ps,req_i1)  
    144105    CALL transfert_request(f_theta_rhodz,req_i1)  
    145106    CALL transfert_request(f_u,req_e1) 
    146      
    147     DO ind=1,ndomain 
    148       CALL swap_dimensions(ind) 
    149       CALL swap_geometry(ind) 
    150        
    151       phis=f_phis(ind) 
    152       ps=f_ps(ind) 
    153       dps=f_dps(ind) 
    154       theta_rhodz=f_theta_rhodz(ind) 
    155       dtheta_rhodz=f_dtheta_rhodz(ind) 
    156       rhodz=f_rhodz(ind) 
    157       p=f_p(ind) 
    158       qu=f_qu(ind) 
    159       u=f_u(ind) 
    160       du=f_du(ind) 
    161       out_u=f_out_u(ind)  
    162 !$OMP PARALLEL DEFAULT(SHARED) 
    163       CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du) 
    164 !$OMP END PARALLEL 
    165     ENDDO 
    166  
    167     IF (mod(it,itau_out)==0 ) THEN 
    168        PRINT *,'CALL write_output_fields' 
    169        CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, & 
    170             f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p) 
    171     END IF 
    172  
    173 !    CALL check_mass_conservation(f_ps,f_dps) 
    174  
    175   END SUBROUTINE caldyn 
    176    
    177    
    178   SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du) 
     107 
     108    SELECT CASE(caldyn_conserv) 
     109    CASE(1) ! energy-conserving 
     110    CASE(2) ! enstrophy-conserving 
     111       DO ind=1,ndomain 
     112          CALL swap_dimensions(ind) 
     113          CALL swap_geometry(ind) 
     114           
     115          phis=f_phis(ind) 
     116          ps=f_ps(ind) 
     117          dps=f_dps(ind) 
     118          theta_rhodz=f_theta_rhodz(ind) 
     119          dtheta_rhodz=f_dtheta_rhodz(ind) 
     120          rhodz=f_rhodz(ind) 
     121          p=f_p(ind) 
     122          qu=f_qu(ind) 
     123          u=f_u(ind) 
     124          du=f_du(ind) 
     125          out_u=f_out_u(ind)  
     126          !$OMP PARALLEL DEFAULT(SHARED) 
     127!          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du) 
     128          CALL compute_pvort(ps, u, p,rhodz,qu) 
     129          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du) 
     130          !$OMP END PARALLEL 
     131       ENDDO 
     132        
     133       IF (mod(it,itau_out)==0 ) THEN 
     134          PRINT *,'CALL write_output_fields' 
     135          CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, & 
     136               f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p) 
     137       END IF 
     138     
     139    CASE DEFAULT 
     140       STOP 
     141    END SELECT 
     142 
     143    !    CALL check_mass_conservation(f_ps,f_dps) 
     144     
     145END SUBROUTINE caldyn 
     146   
     147   
     148SUBROUTINE compute_pvort(ps, u, p,rhodz,qu) 
    179149  USE icosa 
    180150  USE disvert_mod 
    181151  USE exner_mod 
    182152  IMPLICIT NONE 
    183     REAL(rstd),INTENT(IN)  :: phis(iim*jjm) 
    184153    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
    185     REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm) 
    186154    REAL(rstd),INTENT(IN)  :: ps(iim*jjm) 
    187  
    188155    REAL(rstd),INTENT(OUT) :: p(iim*jjm,llm+1) 
    189156    REAL(rstd),INTENT(OUT) :: rhodz(iim*jjm,llm) 
    190157    REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm) 
    191  
    192     REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) 
    193     REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) 
    194     REAL(rstd),INTENT(OUT) :: dps(iim*jjm) 
    195158     
    196159    INTEGER :: i,j,ij,l 
    … …  
    200163 
    201164    REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature 
    202     REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:)   ! Exner function 
    203     REAL(rstd),ALLOCATABLE,SAVE :: pks(:) 
    204 ! Intermediate variable to compute exner function 
    205     REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:) 
    206     REAL(rstd),ALLOCATABLE,SAVE :: beta(:,:) 
    207     
    208     REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential 
    209     REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:)   ! mass flux    
    210     REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux    
    211     REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)    ! mass flux convergence 
    212     REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)       ! vertical velocity       
    213     REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)       ! potential velocity   
    214     REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)   ! bernouilli term 
     165    REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:), pks(:) ! Exner function 
     166    REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:), beta(:,:) 
     167    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)              ! geopotential 
     168    REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:), Ftheta(:,:)  ! mass flux, theta flux 
     169    REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)  ! mass flux convergence 
     170    REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)      ! vertical velocity       
     171    REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)     ! potential velocity   
     172    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function 
    215173     
    216174    LOGICAL,SAVE :: first=.TRUE. 
    … …  
    220178!$OMP MASTER   
    221179!    IF (first) THEN 
    222       ALLOCATE(theta(iim*jjm,llm))    ! potential temperature 
    223       ALLOCATE(pk(iim*jjm,llm))       ! Exner function 
    224       ALLOCATE(pks(iim*jjm)) 
    225       ALLOCATE(alpha(iim*jjm,llm)) 
    226       ALLOCATE(beta(iim*jjm,llm)) 
    227       ALLOCATE(phi(iim*jjm,llm))      ! geopotential 
    228       ALLOCATE(Fe(3*iim*jjm,llm))     ! mass flux    
    229       ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux    
    230       ALLOCATE(convm(iim*jjm,llm))    ! mass flux convergence 
    231       ALLOCATE(w(iim*jjm,llm))        ! vertical velocity       
    232       ALLOCATE(qv(2*iim*jjm,llm))     ! potential velocity   
    233       ALLOCATE(berni(iim*jjm,llm))    ! bernouilli term 
    234 !      first=.FALSE. 
    235 !    ENDIF 
    236 !$OMP END MASTER 
    237 !$OMP BARRIER         
    238  !    du(:,:)=0 
    239 !    theta=1e10 
    240 !    p=1e10 
    241 !    pk=1e10 
    242 !    pks=1e10 
    243 !    alpha=1e10 
    244 !    beta=1e10 
    245 !    phi=1e10 
    246 !    mass=1e10 
    247 !    rhodz=1e10 
    248 !    Fe=1e10 
    249 !    Ftheta=1e10 
    250 !    convm=1e10 
    251 !    w=1e10 
    252 !    qv=1e10 
    253 !    berni=1e10 
    254      
    255  !!! Compute pressure 
    256  
    257 !    PRINT *, 'Computing pressure' 
     180    ALLOCATE(theta(iim*jjm,llm))    ! potential temperature 
     181    ALLOCATE(pk(iim*jjm,llm))       ! Exner function 
     182    ALLOCATE(pks(iim*jjm)) 
     183    ALLOCATE(alpha(iim*jjm,llm)) 
     184    ALLOCATE(beta(iim*jjm,llm)) 
     185    ALLOCATE(phi(iim*jjm,llm))      ! geopotential 
     186    ALLOCATE(Fe(3*iim*jjm,llm))     ! mass flux    
     187    ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux    
     188    ALLOCATE(convm(iim*jjm,llm))    ! mass flux convergence 
     189    ALLOCATE(w(iim*jjm,llm))        ! vertical velocity       
     190    ALLOCATE(qv(2*iim*jjm,llm))     ! potential velocity   
     191    ALLOCATE(berni(iim*jjm,llm))    ! bernouilli term 
     192     
     193!!! Compute pressure 
    258194    DO    l    = 1, llm+1 
    259 !$OMP DO 
    260       DO j=jj_begin-1,jj_end+1 
    261         DO i=ii_begin-1,ii_end+1 
    262           ij=(j-1)*iim+i 
    263           p(ij,l) = ap(l) + bp(l) * ps(ij) 
    264         ENDDO 
    265       ENDDO 
    266     ENDDO 
    267  
     195       !$OMP DO 
     196       DO j=jj_begin-1,jj_end+1 
     197          DO i=ii_begin-1,ii_end+1 
     198             ij=(j-1)*iim+i 
     199             p(ij,l) = ap(l) + bp(l) * ps(ij) 
     200          ENDDO 
     201       ENDDO 
     202    ENDDO 
     203     
    268204!!! Compute mass 
    269    ! PRINT *, 'Computing mass, theta' 
    270205   DO l = 1, llm 
    271206!$OMP DO 
    … …  
    274209            ij=(j-1)*iim+i 
    275210            rhodz(ij,l) = ( p(ij,l) - p(ij,l+1) )/g 
    276             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    277211         ENDDO 
    278212      ENDDO 
    … …  
    309243      ENDDO 
    310244 
    311     DO j=jj_begin,jj_end 
    312       DO i=ii_begin,ii_end 
    313         ij=(j-1)*iim+i 
    314         qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))   
    315         qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))   
    316         qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))   
    317      END DO 
    318   END DO 
    319    
    320 ENDDO 
     245      DO j=jj_begin,jj_end 
     246         DO i=ii_begin,ii_end 
     247            ij=(j-1)*iim+i 
     248            qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))   
     249            qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))   
     250            qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))   
     251         END DO 
     252      END DO 
     253       
     254   ENDDO 
     255 
     256!!$OMP BARRIER 
     257!!$OMP MASTER   
     258    DEALLOCATE(theta)  ! potential temperature 
     259    DEALLOCATE(pk)   ! Exner function 
     260    DEALLOCATE(pks) 
     261    DEALLOCATE(alpha) 
     262    DEALLOCATE(beta) 
     263    DEALLOCATE(phi)    ! geopotential 
     264    DEALLOCATE(Fe)   ! mass flux    
     265    DEALLOCATE(Ftheta) ! theta flux    
     266    DEALLOCATE(convm)    ! mass flux convergence 
     267    DEALLOCATE(w)       ! vertical velocity       
     268    DEALLOCATE(qv)       ! potential velocity   
     269    DEALLOCATE(berni)   ! bernouilli term 
     270!!$OMP END MASTER 
     271!!$OMP BARRIER                                                       
     272  END SUBROUTINE compute_pvort 
     273   
    321274 
    322275!----------- needed next : rhodz, p --------------! 
     276 
     277  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du) 
     278  USE icosa 
     279  USE disvert_mod 
     280  USE exner_mod 
     281  IMPLICIT NONE 
     282    REAL(rstd),INTENT(IN)  :: phis(iim*jjm) 
     283    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
     284    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm) 
     285    REAL(rstd),INTENT(IN)  :: ps(iim*jjm) 
     286    REAL(rstd),INTENT(IN)  :: p(iim*jjm,llm+1) 
     287    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
     288    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm) 
     289 
     290    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) 
     291    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) 
     292    REAL(rstd),INTENT(OUT) :: dps(iim*jjm) 
     293     
     294    INTEGER :: i,j,ij,l 
     295    REAL(rstd) :: ww,uu  
     296    REAL(rstd) :: delta 
     297    REAL(rstd) :: etav,hv, du2 
     298 
     299    REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature 
     300    REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:), pks(:) ! Exner function 
     301    REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:), beta(:,:) 
     302    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)              ! geopotential 
     303    REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:), Ftheta(:,:)  ! mass flux, theta flux 
     304    REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)  ! mass flux convergence 
     305    REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)      ! vertical velocity       
     306    REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)     ! potential velocity   
     307    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function 
     308     
     309    LOGICAL,SAVE :: first=.TRUE. 
     310!$OMP THREADPRIVATE(first) 
     311         
     312!$OMP BARRIER       
     313!$OMP MASTER   
     314!    IF (first) THEN 
     315    ALLOCATE(theta(iim*jjm,llm))    ! potential temperature 
     316    ALLOCATE(pk(iim*jjm,llm))       ! Exner function 
     317    ALLOCATE(pks(iim*jjm)) 
     318    ALLOCATE(alpha(iim*jjm,llm)) 
     319    ALLOCATE(beta(iim*jjm,llm)) 
     320    ALLOCATE(phi(iim*jjm,llm))      ! geopotential 
     321    ALLOCATE(Fe(3*iim*jjm,llm))     ! mass flux    
     322    ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux    
     323    ALLOCATE(convm(iim*jjm,llm))    ! mass flux convergence 
     324    ALLOCATE(w(iim*jjm,llm))        ! vertical velocity       
     325    ALLOCATE(qv(2*iim*jjm,llm))     ! potential velocity   
     326    ALLOCATE(berni(iim*jjm,llm))    ! bernouilli term 
     327 
     328!!! Compute theta 
     329   DO l = 1, llm 
     330!$OMP DO 
     331      DO j=jj_begin-1,jj_end+1 
     332         DO i=ii_begin-1,ii_end+1 
     333            ij=(j-1)*iim+i 
     334            theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
     335         ENDDO 
     336      ENDDO 
     337   ENDDO 
    323338 
    324339  DO l = 1, llm 
    … …  
    577592!!$OMP BARRIER                                                       
    578593  END SUBROUTINE compute_caldyn 
     594 
     595!-------------------------------- Diagnostics ---------------------------- 
     596 
     597  SUBROUTINE check_mass_conservation(f_ps,f_dps) 
     598  USE icosa 
     599  IMPLICIT NONE 
     600    TYPE(t_field),POINTER :: f_ps(:) 
     601    TYPE(t_field),POINTER :: f_dps(:) 
     602    REAL(rstd),POINTER :: ps(:) 
     603    REAL(rstd),POINTER :: dps(:) 
     604    REAL(rstd) :: mass_tot,dmass_tot 
     605    INTEGER :: ind,i,j,ij 
     606     
     607    mass_tot=0 
     608    dmass_tot=0 
     609     
     610    CALL transfert_request(f_dps,req_i1) 
     611    CALL transfert_request(f_ps,req_i1) 
     612 
     613    DO ind=1,ndomain 
     614      CALL swap_dimensions(ind) 
     615      CALL swap_geometry(ind) 
     616 
     617      ps=f_ps(ind) 
     618      dps=f_dps(ind) 
     619 
     620      DO j=jj_begin,jj_end 
     621        DO i=ii_begin,ii_end 
     622          ij=(j-1)*iim+i 
     623          IF (domain(ind)%own(i,j)) THEN 
     624            mass_tot=mass_tot+ps(ij)*Ai(ij)/g 
     625            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g 
     626          ENDIF 
     627        ENDDO 
     628      ENDDO 
     629    
     630    ENDDO 
     631    PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot         
     632 
     633  END SUBROUTINE check_mass_conservation   
    579634   
    580635  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, & 
Note: See TracChangeset for help on using the changeset viewer.