Changeset 123


Ignore:
Timestamp:
02/01/13 00:52:47 (11 years ago)
Author:
dubos
Message:

Some reorganization of caldyn.gcm

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r122 r123  
    179179    REAL(rstd) :: etav,hv, du2 
    180180 
    181 !   REAL(rstd) :: theta(iim*jjm,llm)  ! potential temperature 
    182 !   REAL(rstd) :: p(iim*jjm,llm+1)  ! pression 
    183 !   REAL(rstd) :: pk(iim*jjm,llm)   ! Exner function 
    184 !   REAL(rstd) :: pks(iim*jjm) 
    185 !! Intermediate variable to compute exner function 
    186 !   REAL(rstd) :: alpha(iim*jjm,llm) 
    187 !   REAL(rstd) :: beta(iim*jjm,llm) 
    188 !!    
    189 !   REAL(rstd) :: phi(iim*jjm,llm)    ! geopotential 
    190 !   REAL(rstd) :: mass(iim*jjm,llm)   ! mass     
    191 !   REAL(rstd) :: rhodz(iim*jjm,llm)   ! mass density    
    192 !   REAL(rstd) :: Fe(3*iim*jjm,llm)   ! mass flux    
    193 !   REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux    
    194 !   REAL(rstd) :: convm(iim*jjm,llm)    ! mass flux convergence 
    195 !   REAL(rstd) :: w(iim*jjm,llm)       ! vertical velocity       
    196 !   REAL(rstd) :: qv(2*iim*jjm,llm)       ! potential velocity   
    197 !   REAL(rstd) :: berni(iim*jjm,llm)   ! bernouilli term 
    198      
    199181    REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature 
    200182    REAL(rstd),ALLOCATABLE,SAVE :: p(:,:)  ! pression 
    … …  
    206188 !    
    207189    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential 
    208     REAL(rstd),ALLOCATABLE,SAVE :: mass(:,:)   ! mass     
    209190    REAL(rstd),ALLOCATABLE,SAVE :: rhodz(:,:)   ! mass density    
    210191    REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:)   ! mass flux    
    … …  
    228209      ALLOCATE(beta(iim*jjm,llm)) 
    229210      ALLOCATE(phi(iim*jjm,llm))    ! geopotential 
    230       ALLOCATE(mass(iim*jjm,llm))   ! mass     
    231211      ALLOCATE(rhodz(iim*jjm,llm))   ! mass density    
    232212      ALLOCATE(Fe(3*iim*jjm,llm))   ! mass flux    
    … …  
    269249      ENDDO 
    270250    ENDDO 
    271        
    272  !!! Compute Exner function 
    273 !    PRINT *, 'Computing Exner' 
    274     CALL compute_exner(ps,p,pks,pk,1)  
    275251 
    276252!!! Compute mass 
    277    ! PRINT *, 'Computing mass' 
     253   ! PRINT *, 'Computing mass, theta' 
    278254   DO l = 1, llm 
    279255!$OMP DO 
    280      DO j=jj_begin-1,jj_end+1 
    281        DO i=ii_begin-1,ii_end+1 
    282          ij=(j-1)*iim+i 
    283          mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g 
    284          rhodz(ij,l) = mass(ij,l) / Ai(ij) 
    285        ENDDO 
    286      ENDDO 
     256      DO j=jj_begin-1,jj_end+1 
     257         DO i=ii_begin-1,ii_end+1 
     258            ij=(j-1)*iim+i 
     259            rhodz(ij,l) = ( p(ij,l) - p(ij,l+1) )/g 
     260            theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
     261         ENDDO 
     262      ENDDO 
    287263   ENDDO 
    288  
    289 !! compute theta 
    290    ! PRINT *, 'Computing theta' 
    291     DO    l    = 1, llm 
    292 !$OMP DO 
    293       DO j=jj_begin-1,jj_end+1 
    294         DO i=ii_begin-1,ii_end+1 
    295           ij=(j-1)*iim+i 
    296           theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    297         ENDDO 
    298       ENDDO 
    299     ENDDO 
    300  
    301  
    302 !!! Compute geopotential 
    303  
    304   ! for first layer 
    305 !$OMP DO 
    306    DO j=jj_begin-1,jj_end+1 
    307      DO i=ii_begin-1,ii_end+1 
    308        ij=(j-1)*iim+i 
    309        phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) ) 
    310      ENDDO 
    311    ENDDO 
    312            
    313   ! for other layers 
    314    DO l = 2, llm 
    315 !$OMP DO 
    316      DO j=jj_begin-1,jj_end+1 
    317        DO i=ii_begin-1,ii_end+1 
    318          ij=(j-1)*iim+i 
    319          phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  &  
    320                                        * (  pk(ij,l-1) -  pk(ij,l)    ) 
    321        ENDDO 
    322      ENDDO 
    323    ENDDO 
    324  
    325     
    326 !!!  Compute mass flux 
    327  
    328   DO l = 1, llm 
    329 !$OMP DO 
    330     DO j=jj_begin-1,jj_end+1 
    331       DO i=ii_begin-1,ii_end+1 
    332         ij=(j-1)*iim+i 
    333         Fe(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 
    334         Fe(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 
    335         Fe(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 
    336       ENDDO 
    337     ENDDO 
    338   ENDDO 
    339  
    340 !!! fisrt composante dtheta 
    341  
    342  ! Flux on the edge  
    343   DO l = 1, llm 
    344 !$OMP DO 
    345     DO j=jj_begin-1,jj_end+1 
    346       DO i=ii_begin-1,ii_end+1 
    347         ij=(j-1)*iim+i   
    348         Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*Fe(ij+u_right,l) 
    349         Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*Fe(ij+u_lup,l) 
    350         Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*Fe(ij+u_ldown,l) 
    351       ENDDO 
    352     ENDDO 
    353   ENDDO 
    354   
    355   
    356  ! compute divergence  
    357   DO l = 1, llm 
    358 !$OMP DO 
    359     DO j=jj_begin,jj_end 
    360       DO i=ii_begin,ii_end 
    361         ij=(j-1)*iim+i   
    362 ! signe ? attention d (rho theta dz) 
    363         ! dtheta_rhodz =  -div(flux.theta) 
    364         dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  & 
    365                                  ne(ij,rup)*Ftheta(ij+u_rup,l)        +  &   
    366                                  ne(ij,lup)*Ftheta(ij+u_lup,l)        +  &   
    367                                  ne(ij,left)*Ftheta(ij+u_left,l)      +  &   
    368                                  ne(ij,ldown)*Ftheta(ij+u_ldown,l)    +  &   
    369                                  ne(ij,rdown)*Ftheta(ij+u_rdown,l)) 
    370       ENDDO 
    371     ENDDO 
    372   ENDDO 
    373  
    374  
    375  
    376 !!! mass flux convergence computation   
    377  
    378  ! horizontal convergence 
    379   DO l = 1, llm 
    380 !$OMP DO 
    381     DO j=jj_begin,jj_end 
    382       DO i=ii_begin,ii_end 
    383         ij=(j-1)*iim+i 
    384         ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 
    385         convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l)   +  & 
    386                                 ne(ij,rup)*Fe(ij+u_rup,l)       +  &   
    387                                 ne(ij,lup)*Fe(ij+u_lup,l)       +  &   
    388                                 ne(ij,left)*Fe(ij+u_left,l)     +  &   
    389                                 ne(ij,ldown)*Fe(ij+u_ldown,l)   +  &   
    390                                 ne(ij,rdown)*Fe(ij+u_rdown,l))     
    391        ENDDO 
    392      ENDDO 
    393    ENDDO 
    394   
    395     
    396  ! vertical integration from up to down 
    397   DO  l = llm-1, 1, -1 
    398 !$OMP DO 
    399     DO j=jj_begin,jj_end 
    400       DO i=ii_begin,ii_end 
    401         ij=(j-1)*iim+i 
    402         convm(ij,l) = convm(ij,l) + convm(ij,l+1) 
    403       ENDDO 
    404     ENDDO 
    405   ENDDO 
    406  
    407    
    408 !!! Compute dps 
    409 !$OMP DO 
    410   DO j=jj_begin,jj_end 
    411     DO i=ii_begin,ii_end 
    412       ij=(j-1)*iim+i 
    413       ! dps/dt = -int(div flux)dz 
    414       dps(ij)=-convm(ij,1) * g  
    415     ENDDO 
    416   ENDDO 
    417  
    418  
    419 !!! Compute vertical velocity 
    420   DO l = 1,llm-1 
    421 !$OMP DO 
    422     DO j=jj_begin,jj_end 
    423       DO i=ii_begin,ii_end 
    424         ij=(j-1)*iim+i 
    425         ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 
    426         ! => w>0 for upward transport 
    427         w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ) 
    428       ENDDO 
    429     ENDDO 
    430   ENDDO 
    431  
    432 !$OMP DO 
    433   ! vertical mass flux at the surface = 0 
    434   DO j=jj_begin,jj_end 
    435     DO i=ii_begin,ii_end 
    436       ij=(j-1)*iim+i 
    437       w(ij,1)  = 0. 
    438     ENDDO 
    439   ENDDO 
    440  
    441264 
    442265!!! Compute shallow-water potential vorticity 
    … …  
    471294    ENDDO 
    472295 
     296 
     297  DO l = 1, llm 
     298!!!  Compute mass and theta fluxes 
     299!$OMP DO 
     300    DO j=jj_begin-1,jj_end+1 
     301      DO i=ii_begin-1,ii_end+1 
     302        ij=(j-1)*iim+i 
     303        Fe(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 
     304        Fe(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 
     305        Fe(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 
     306 
     307        Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*Fe(ij+u_right,l) 
     308        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*Fe(ij+u_lup,l) 
     309        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*Fe(ij+u_ldown,l) 
     310      ENDDO 
     311    ENDDO 
     312 
     313!!! compute horizontal divergence of fluxes 
     314!$OMP DO 
     315    DO j=jj_begin,jj_end 
     316      DO i=ii_begin,ii_end 
     317        ij=(j-1)*iim+i   
     318        ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 
     319        convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l)   +  & 
     320                                ne(ij,rup)*Fe(ij+u_rup,l)       +  &   
     321                                ne(ij,lup)*Fe(ij+u_lup,l)       +  &   
     322                                ne(ij,left)*Fe(ij+u_left,l)     +  &   
     323                                ne(ij,ldown)*Fe(ij+u_ldown,l)   +  &   
     324                                ne(ij,rdown)*Fe(ij+u_rdown,l))     
     325 
     326        ! signe ? attention d (rho theta dz) 
     327        ! dtheta_rhodz =  -div(flux.theta) 
     328        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  & 
     329                                 ne(ij,rup)*Ftheta(ij+u_rup,l)        +  &   
     330                                 ne(ij,lup)*Ftheta(ij+u_lup,l)        +  &   
     331                                 ne(ij,left)*Ftheta(ij+u_left,l)      +  &   
     332                                 ne(ij,ldown)*Ftheta(ij+u_ldown,l)    +  &   
     333                                 ne(ij,rdown)*Ftheta(ij+u_rdown,l)) 
     334      ENDDO 
     335    ENDDO 
     336  ENDDO 
     337    
     338!!! cumulate mass flux convergence from top to bottom 
     339  DO  l = llm-1, 1, -1 
     340!$OMP DO 
     341    DO j=jj_begin,jj_end 
     342      DO i=ii_begin,ii_end 
     343        ij=(j-1)*iim+i 
     344        convm(ij,l) = convm(ij,l) + convm(ij,l+1) 
     345      ENDDO 
     346    ENDDO 
     347  ENDDO 
     348   
     349!!! Compute vertical mass flux 
     350  DO l = 1,llm-1 
     351!$OMP DO 
     352    DO j=jj_begin,jj_end 
     353      DO i=ii_begin,ii_end 
     354        ij=(j-1)*iim+i 
     355        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 
     356        ! => w>0 for upward transport 
     357        w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ) 
     358      ENDDO 
     359    ENDDO 
     360  ENDDO 
     361 
     362! compute dps, vertical mass flux at the surface = 0 
     363!$OMP DO 
     364  DO j=jj_begin,jj_end 
     365    DO i=ii_begin,ii_end 
     366      ij=(j-1)*iim+i 
     367      w(ij,1)  = 0. 
     368      ! dps/dt = -int(div flux)dz 
     369      dps(ij)=-convm(ij,1) * g  
     370    ENDDO 
     371  ENDDO 
     372 
    473373!!! Compute potential vorticity (Coriolis) contribution to du 
    474374  DO l=1,llm 
    … …  
    521421  ENDDO 
    522422   
    523  
     423!!! Compute Exner function 
     424!    PRINT *, 'Computing Exner' 
     425    CALL compute_exner(ps,p,pks,pk,1)  
     426 
     427!!! Compute geopotential 
     428 
     429  ! for first layer 
     430!$OMP DO 
     431   DO j=jj_begin-1,jj_end+1 
     432     DO i=ii_begin-1,ii_end+1 
     433       ij=(j-1)*iim+i 
     434       phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) ) 
     435     ENDDO 
     436   ENDDO 
     437           
     438  ! for other layers 
     439   DO l = 2, llm 
     440!$OMP DO 
     441     DO j=jj_begin-1,jj_end+1 
     442       DO i=ii_begin-1,ii_end+1 
     443         ij=(j-1)*iim+i 
     444         phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  &  
     445                                       * (  pk(ij,l-1) -  pk(ij,l)    ) 
     446       ENDDO 
     447     ENDDO 
     448   ENDDO 
     449 
     450    
    524451!!! Compute bernouilli term = Kinetic Energy + geopotential 
    525452  DO l=1,llm 
    … …  
    542469   
    543470  
    544 !!! second contribution to du (gradients of Bernoulli and Exner functions) 
    545   DO l=1,llm 
    546 !$OMP DO 
    547     DO j=jj_begin,jj_end 
    548       DO i=ii_begin,ii_end 
    549         ij=(j-1)*iim+i 
    550          
    551          du(ij+u_right,l)= du(ij+u_right,l)+ 1/de(ij+u_right) * (                              & 
    552                            0.5*(theta(ij,l)+theta(ij+t_right,l))                             & 
    553                               *( ne(ij,right)*pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) & 
    554                          + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) ) 
    555           
    556          du(ij+u_lup,l)= du(ij+u_lup,l)+ 1/de(ij+u_lup) * (                                  & 
    557                            0.5*(theta(ij,l)+theta(ij+t_lup,l))                             & 
    558                               *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l))    & 
    559                          + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) ) 
    560                           
    561          du(ij+u_ldown,l)= du(ij+u_ldown,l)+ 1/de(ij+u_ldown) * (                                & 
    562                            0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               & 
    563                               *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l))    & 
    564                          + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) ) 
    565       ENDDO 
    566     ENDDO 
    567   ENDDO 
    568   
    569 !!! save second contribution to du for debugging output 
     471!!! gradients of Bernoulli and Exner functions 
    570472  DO l=1,llm 
    571473!$OMP DO 
    … …  
    579481                        + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) ) 
    580482         
     483        du(ij+u_right,l) = du(ij+u_right,l) + out_u(ij+u_right,l) 
     484         
    581485        out_u(ij+u_lup,l)=  1/de(ij+u_lup) * (                                  & 
    582486                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                             & 
    583487                             *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l))    & 
    584488                        + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) ) 
    585                          
     489         
     490        du(ij+u_lup,l) = du(ij+u_lup,l) + out_u(ij+u_lup,l) 
     491                 
    586492        out_u(ij+u_ldown,l)= 1/de(ij+u_ldown) * (                                & 
    587493                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               & 
    588494                             *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l))    & 
    589495                        + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) ) 
     496         
     497        du(ij+u_ldown,l) = du(ij+u_ldown,l) + out_u(ij+u_ldown,l) 
     498 
    590499      ENDDO 
    591500    ENDDO 
    592501  ENDDO 
    593  
    594 !!! contributions due to vertical advection 
    595502  
    596  ! Contribution to dtheta 
     503!!! contributions of vertical advection to du, dtheta 
     504  
    597505  DO l=1,llm-1 
    598506!$OMP DO 
    599507    DO j=jj_begin,jj_end 
    600508      DO i=ii_begin,ii_end 
     509        ij=(j-1)*iim+i 
    601510         ! ww>0 <=> upward transport 
    602         ij=(j-1)*iim+i 
     511 
    603512        ww            =   0.5 * w(ij,l+1) * (theta(ij,l) + theta(ij,l+1) ) 
    604513        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww 
    605514        dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww 
    606       ENDDO 
    607     ENDDO 
    608   ENDDO        
    609  
    610  
    611  ! Contribution to du 
    612   DO l=1,llm-1 
    613 !$OMP DO 
    614     DO j=jj_begin,jj_end 
    615       DO i=ii_begin,ii_end 
    616         ij=(j-1)*iim+i 
     515 
    617516        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_right,l+1)) 
    618517        uu  = u(ij+u_right,l+1) - u(ij+u_right,l)   
    … …  
    643542    DEALLOCATE(beta) 
    644543    DEALLOCATE(phi)    ! geopotential 
    645     DEALLOCATE(mass)   ! mass     
     544!    DEALLOCATE(mass)   ! mass     
    646545    DEALLOCATE(rhodz)   ! mass density    
    647546    DEALLOCATE(Fe)   ! mass flux    
Note: See TracChangeset for help on using the changeset viewer.