Changeset 361 for codes/icosagcm/trunk


Ignore:
Timestamp:
09/09/15 15:26:13 (9 years ago)
Author:
dubos
Message:

ARK2.3 HEVI time stepping - tested with DCMIP4 only

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

Legend:

Unmodified
Added
Removed

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

    r360 r361  
    22  USE icosa 
    33  USE transfert_mod 
     4  USE caldyn_kernels_mod 
    45  IMPLICIT NONE 
    56  PRIVATE 
    6  
    7   INTEGER, PARAMETER :: energy=1, enstrophy=2 
    8   TYPE(t_field),POINTER :: f_out_u(:), f_qu(:), f_qv(:) 
    9   REAL(rstd),SAVE,POINTER :: out_u(:,:), p(:,:), qu(:,:) 
    10   !$OMP THREADPRIVATE(out_u, p, qu) 
    11  
    12 ! temporary shared variable for caldyn 
    13   TYPE(t_field),POINTER :: f_pk(:) 
    14   TYPE(t_field),POINTER :: f_wwuu(:) 
    15   TYPE(t_field),POINTER :: f_planetvel(:) 
    16     
    17   INTEGER :: caldyn_conserv 
    18  !$OMP THREADPRIVATE(caldyn_conserv)  
    19   
    20   TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu 
    217 
    228  PUBLIC init_caldyn, caldyn_BC, caldyn, req_ps, req_mass 
    … …  
    163149!$OMP THREADPRIVATE(first) 
    164150     
    165     ! MPI messages need to be sent at first call to caldyn 
    166     ! This is needed only once : the next ones will be sent by timeloop 
    167151    IF (first) THEN  
    168152      first=.FALSE. 
    … …  
    175159      CALL init_message(f_u,req_e1_vect,req_u) 
    176160      CALL init_message(f_qu,req_e1_scal,req_qu) 
     161      ! Overlapping com/compute (deactivated) : MPI messages need to be sent at first call to caldyn 
     162      ! This is needed only once : the next ones will be sent by timeloop 
    177163!      IF(caldyn_eta==eta_mass) THEN 
    178164!         CALL send_message(f_ps,req_ps)  
    … …  
    186172    CALL trace_start("caldyn") 
    187173 
    188       IF(caldyn_eta==eta_mass) THEN 
    189          CALL send_message(f_ps,req_ps)  
    190       ELSE 
    191          CALL send_message(f_mass,req_mass)  
    192       END IF 
    193  
    194     CALL send_message(f_theta_rhodz,req_theta_rhodz)  
    195     CALL send_message(f_u,req_u) 
     174    IF(caldyn_eta==eta_mass) THEN 
     175       CALL send_message(f_ps,req_ps) ! COM00 
     176    ELSE 
     177       CALL send_message(f_mass,req_mass) ! COM00 
     178    END IF 
     179     
     180    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01 
     181    CALL send_message(f_u,req_u) ! COM02 
    196182 
    197183    SELECT CASE(caldyn_conserv) 
    … …  
    208194          qu=f_qu(ind) 
    209195          qv=f_qv(ind) 
    210           CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv) 
     196          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv) ! COM00 COM01 COM02 
    211197       ENDDO 
    212198!       CALL checksum(f_mass) 
    213199!       CALL checksum(f_theta) 
    214200 
    215        CALL send_message(f_qu,req_qu) 
     201       CALL send_message(f_qu,req_qu) ! COM03 wait_message in caldyn_horiz 
    216202!       CALL wait_message(req_qu) 
    217203 
    … …  
    287273END SUBROUTINE caldyn 
    288274 
    289 SUBROUTINE compute_planetvel(planetvel) 
    290   USE wind_mod 
    291   REAL(rstd),INTENT(OUT)  :: planetvel(iim*3*jjm) 
    292   REAL(rstd) :: ulon(iim*3*jjm) 
    293   REAL(rstd) :: ulat(iim*3*jjm) 
    294   REAL(rstd) :: lon,lat 
    295   INTEGER :: ij 
    296   DO ij=ij_begin_ext,ij_end_ext 
    297      ulon(ij+u_right)=radius*omega*cos(lat_e(ij+u_right)) 
    298      ulat(ij+u_right)=0 
    299  
    300      ulon(ij+u_lup)=radius*omega*cos(lat_e(ij+u_lup)) 
    301      ulat(ij+u_lup)=0 
    302  
    303      ulon(ij+u_ldown)=radius*omega*cos(lat_e(ij+u_ldown)) 
    304      ulat(ij+u_ldown)=0 
    305   END DO 
    306   CALL compute_wind2D_perp_from_lonlat_compound(ulon, ulat, planetvel) 
    307 END SUBROUTINE compute_planetvel 
    308      
    309 SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv) 
    310   USE icosa 
    311   USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass 
    312   USE exner_mod 
    313   USE trace 
    314   USE omp_para 
    315   IMPLICIT NONE 
    316   REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
    317   REAL(rstd),INTENT(IN)  :: ps(iim*jjm) 
    318   REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm) 
    319   REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm) 
    320   REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm) 
    321   REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm) 
    322   REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm) 
    323    
    324   INTEGER :: i,j,ij,l 
    325   REAL(rstd) :: etav,hv, m 
    326 !  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity   
    327    
    328   CALL trace_start("compute_pvort")   
    329  
    330   IF(caldyn_eta==eta_mass) THEN 
    331      CALL wait_message(req_ps)   
    332   ELSE 
    333      CALL wait_message(req_mass) 
    334   END IF 
    335   CALL wait_message(req_theta_rhodz)  
    336  
    337   IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta 
    338      DO l = ll_begin,ll_end 
    339         CALL test_message(req_u)  
    340 !DIR$ SIMD 
    341         DO ij=ij_begin_ext,ij_end_ext 
    342            m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 
    343            rhodz(ij,l) = m 
    344            theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    345         ENDDO 
    346      ENDDO 
    347   ELSE ! Compute only theta 
    348      DO l = ll_begin,ll_end 
    349         CALL test_message(req_u)  
    350 !DIR$ SIMD 
    351        DO ij=ij_begin_ext,ij_end_ext 
    352          theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    353        ENDDO 
    354      ENDDO 
    355   END IF 
    356  
    357   CALL wait_message(req_u)    
    358    
    359 !!! Compute shallow-water potential vorticity 
    360   DO l = ll_begin,ll_end 
    361 !DIR$ SIMD 
    362      DO ij=ij_begin_ext,ij_end_ext 
    363           etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         & 
    364                                 + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  & 
    365                                 - ne_lup        * u(ij+u_lup,l)        * de(ij+u_lup) )                                
    366  
    367           hv =  Riv2(ij,vup)          * rhodz(ij,l)            & 
    368               + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     & 
    369               + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 
    370       
    371           qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv 
    372            
    373           etav = 1./Av(ij+z_down)*(  ne_ldown         * u(ij+u_ldown,l)          * de(ij+u_ldown)          & 
    374                                    + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  & 
    375                                    - ne_rdown         * u(ij+u_rdown,l)          * de(ij+u_rdown) ) 
    376         
    377           hv =  Riv2(ij,vdown)        * rhodz(ij,l)          & 
    378               + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  & 
    379               + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 
    380                         
    381           qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 
    382  
    383       ENDDO 
    384  
    385 !DIR$ SIMD 
    386       DO ij=ij_begin,ij_end 
    387          qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))   
    388          qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))   
    389          qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))   
    390       END DO 
    391        
    392    ENDDO 
    393  
    394    CALL trace_end("compute_pvort") 
    395                                                         
    396   END SUBROUTINE compute_pvort 
    397    
    398   SUBROUTINE compute_geopot(ps,rhodz,theta, pk,geopot)  
    399   USE icosa 
    400   USE disvert_mod 
    401   USE exner_mod 
    402   USE trace 
    403   USE omp_para 
    404   IMPLICIT NONE 
    405     REAL(rstd),INTENT(INOUT) :: ps(iim*jjm) 
    406     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm) 
    407     REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)    ! potential temperature 
    408     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)       ! Exner function 
    409     REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential 
    410      
    411     INTEGER :: i,j,ij,l 
    412     REAL(rstd) :: p_ik, exner_ik 
    413     INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext 
    414  
    415  
    416     CALL trace_start("compute_geopot") 
    417      
    418     CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext) 
    419     ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1 
    420     ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1 
    421  
    422     IF(caldyn_eta==eta_mass) THEN 
    423  
    424 !!! Compute exner function and geopotential 
    425          DO l = 1,llm 
    426           !DIR$ SIMD 
    427             DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    428                p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 
    429                !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j)) 
    430                exner_ik = cpp * (p_ik/preff) ** kappa 
    431                pk(ij,l) = exner_ik 
    432              ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
    433                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
    434           ENDDO 
    435          ENDDO 
    436 !       ENDIF 
    437     ELSE  
    438        ! We are using a Lagrangian vertical coordinate 
    439        ! Pressure must be computed first top-down (temporarily stored in pk) 
    440        ! Then Exner pressure and geopotential are computed bottom-up 
    441        ! Notice that the computation below should work also when caldyn_eta=eta_mass           
    442  
    443        IF(boussinesq) THEN ! compute only geopotential : pressure pk will be computed in compute_caldyn_horiz 
    444           ! specific volume 1 = dphi/g/rhodz 
    445 !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency 
    446          DO l = 1,llm 
    447            !DIR$ SIMD 
    448            DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    449               geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l) 
    450            ENDDO 
    451          ENDDO 
    452        ELSE ! non-Boussinesq, compute geopotential and Exner pressure 
    453           ! uppermost layer 
    454  
    455          !DIR$ SIMD 
    456          DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    457             pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 
    458          END DO 
    459          ! other layers 
    460          DO l = llm-1, 1, -1 
    461             !DIR$ SIMD 
    462             DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    463                pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1)) 
    464             END DO 
    465          END DO 
    466         ! surface pressure (for diagnostics) 
    467          DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    468             ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 
    469          END DO 
    470  
    471          ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
    472          DO l = 1,llm 
    473            !DIR$ SIMD 
    474             DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    475                p_ik = pk(ij,l) 
    476                exner_ik = cpp * (p_ik/preff) ** kappa 
    477                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik  
    478                pk(ij,l) = exner_ik 
    479             ENDDO 
    480           ENDDO 
    481        END IF 
    482  
    483     END IF 
    484  
    485 !ym flush geopot 
    486 !$OMP BARRIER 
    487  
    488   CALL trace_end("compute_geopot") 
    489  
    490   END SUBROUTINE compute_geopot 
    491  
    492   SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du) 
    493   USE icosa 
    494   USE disvert_mod 
    495   USE exner_mod 
    496   USE trace 
    497   USE omp_para 
    498   IMPLICIT NONE 
    499     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)    ! prognostic "velocity" 
    500     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    501     REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm) 
    502     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature 
    503     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function 
    504     REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential 
    505  
    506     REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s 
    507     REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence 
    508     REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) 
    509     REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) 
    510  
    511     REAL(rstd) :: cor_NT(iim*jjm,llm)  ! NT coriolis force u.(du/dPhi) 
    512     REAL(rstd) :: urel(3*iim*jjm,llm)  ! relative velocity 
    513     REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux 
    514     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
    515      
    516     INTEGER :: i,j,ij,l 
    517     REAL(rstd) :: ww,uu 
    518  
    519     CALL trace_start("compute_caldyn_horiz") 
    520  
    521 !    CALL wait_message(req_theta_rhodz)  
    522  
    523   DO l = ll_begin, ll_end 
    524 !!!  Compute mass and theta fluxes 
    525     IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    526 !DIR$ SIMD 
    527     DO ij=ij_begin_ext,ij_end_ext          
    528         hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 
    529         hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 
    530         hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 
    531  
    532         Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l) 
    533         Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l) 
    534         Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l) 
    535     ENDDO 
    536  
    537 !!! compute horizontal divergence of fluxes 
    538 !DIR$ SIMD 
    539     DO ij=ij_begin,ij_end          
    540         ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 
    541         convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  & 
    542                                 ne_rup*hflux(ij+u_rup,l)       +  &   
    543                                 ne_lup*hflux(ij+u_lup,l)       +  &   
    544                                 ne_left*hflux(ij+u_left,l)     +  &   
    545                                 ne_ldown*hflux(ij+u_ldown,l)   +  &   
    546                                 ne_rdown*hflux(ij+u_rdown,l))     
    547  
    548         ! signe ? attention d (rho theta dz) 
    549         ! dtheta_rhodz =  -div(flux.theta) 
    550         dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  & 
    551                                  ne_rup*Ftheta(ij+u_rup,l)        +  &   
    552                                  ne_lup*Ftheta(ij+u_lup,l)        +  &   
    553                                  ne_left*Ftheta(ij+u_left,l)      +  &   
    554                                  ne_ldown*Ftheta(ij+u_ldown,l)    +  &   
    555                                  ne_rdown*Ftheta(ij+u_rdown,l)) 
    556     ENDDO 
    557  
    558  END DO 
    559  
    560 !!! Compute potential vorticity (Coriolis) contribution to du 
    561   
    562     SELECT CASE(caldyn_conserv) 
    563     CASE(energy) ! energy-conserving TRiSK 
    564  
    565        CALL wait_message(req_qu) 
    566  
    567         DO l=ll_begin,ll_end 
    568 !DIR$ SIMD 
    569           DO ij=ij_begin,ij_end          
    570      
    571                 uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             & 
    572                      wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             & 
    573                      wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           & 
    574                      wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         & 
    575                      wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         &  
    576                      wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         & 
    577                      wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         & 
    578                      wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         & 
    579                      wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             & 
    580                      wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l)) 
    581                 du(ij+u_right,l) = .5*uu/de(ij+u_right) 
    582                  
    583                 uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        & 
    584                      wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      & 
    585                      wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      & 
    586                      wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      & 
    587                      wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          &  
    588                      wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          & 
    589                      wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              & 
    590                      wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              & 
    591                      wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            & 
    592                      wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l)) 
    593                 du(ij+u_lup,l) = .5*uu/de(ij+u_lup) 
    594  
    595                  
    596                 uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      & 
    597                      wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      & 
    598                      wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          & 
    599                      wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          & 
    600                      wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        &  
    601                      wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          & 
    602                      wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        & 
    603                      wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      & 
    604                      wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      & 
    605                      wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l)) 
    606                 du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown) 
    607                  
    608           ENDDO 
    609        ENDDO 
    610         
    611     CASE(enstrophy) ! enstrophy-conserving TRiSK 
    612    
    613         DO l=ll_begin,ll_end 
    614 !DIR$ SIMD 
    615           DO ij=ij_begin,ij_end          
    616  
    617                 uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           & 
    618                      wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           & 
    619                      wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          & 
    620                      wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         & 
    621                      wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         &  
    622                      wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 & 
    623                      wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 & 
    624                      wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 & 
    625                      wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   & 
    626                      wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l) 
    627                 du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right) 
    628                  
    629                  
    630                 uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        & 
    631                      wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       & 
    632                      wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       & 
    633                      wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       & 
    634                      wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         &  
    635                      wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 & 
    636                      wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   & 
    637                      wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   & 
    638                      wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  & 
    639                      wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l) 
    640                 du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup) 
    641  
    642                 uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      & 
    643                      wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      & 
    644                      wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        & 
    645                      wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        & 
    646                      wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       &  
    647                      wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                & 
    648                      wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               & 
    649                      wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              & 
    650                      wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              & 
    651                      wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l) 
    652                 du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown) 
    653       
    654           ENDDO 
    655        ENDDO 
    656         
    657     CASE DEFAULT 
    658        STOP 
    659     END SELECT 
    660    
    661 !!! Compute bernouilli term = Kinetic Energy + geopotential 
    662     IF(boussinesq) THEN 
    663        ! first use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure)  
    664        ! uppermost layer 
    665        !DIR$ SIMD 
    666        DO ij=ij_begin_ext,ij_end_ext          
    667           pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm) 
    668        END DO 
    669        ! other layers 
    670        DO l = llm-1, 1, -1 
    671 !          !$OMP DO SCHEDULE(STATIC)  
    672           !DIR$ SIMD 
    673           DO ij=ij_begin_ext,ij_end_ext          
    674              pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1)) 
    675           END DO 
    676        END DO 
    677        ! surface pressure (for diagnostics) FIXME 
    678        ! DO ij=ij_begin_ext,ij_end_ext          
    679        !   ps(ij) = pk(ij,1) + (.5*g)*theta(ij,1)*rhodz(ij,1) 
    680        ! END DO 
    681        ! now pk contains the Lagrange multiplier (pressure) 
    682  
    683        DO l=ll_begin,ll_end 
    684 !DIR$ SIMD 
    685           DO ij=ij_begin,ij_end          
    686                  
    687                 berni(ij,l) = pk(ij,l) + & 
    688                        1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    & 
    689                                      le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
    690                                      le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
    691                                      le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       & 
    692                                      le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
    693                                      le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    694                 ! from now on pk contains the vertically-averaged geopotential 
    695                 pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 
    696           ENDDO 
    697        ENDDO 
    698  
    699     ELSE ! compressible 
    700  
    701        DO l=ll_begin,ll_end 
    702 !DIR$ SIMD 
    703             DO ij=ij_begin,ij_end          
    704                  
    705                 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
    706                      + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    & 
    707                                      le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
    708                                      le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
    709                                      le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       & 
    710                                      le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
    711                                      le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    712           ENDDO 
    713        ENDDO 
    714  
    715     END IF ! Boussinesq/compressible 
    716  
    717 !!! Add gradients of Bernoulli and Exner functions to du 
    718     DO l=ll_begin,ll_end 
    719 !DIR$ SIMD 
    720        DO ij=ij_begin,ij_end          
    721          
    722              du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             & 
    723                                0.5*(theta(ij,l)+theta(ij+t_right,l))                & 
    724                                   *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    & 
    725                                    + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) ) 
    726          
    727          
    728              du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  & 
    729                                0.5*(theta(ij,l)+theta(ij+t_lup,l))                  & 
    730                                   *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       & 
    731                                    + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) ) 
    732                  
    733              du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             & 
    734                                0.5*(theta(ij,l)+theta(ij+t_ldown,l))                & 
    735                                   *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     & 
    736                                    + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) ) 
    737  
    738        ENDDO 
    739     ENDDO 
    740   
    741     CALL trace_end("compute_caldyn_horiz") 
    742  
    743 END SUBROUTINE compute_caldyn_horiz 
    744  
    745 SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du) 
    746   USE icosa 
    747   USE disvert_mod 
    748   USE exner_mod 
    749   USE trace 
    750   USE omp_para 
    751   IMPLICIT NONE 
    752     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
    753     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm) 
    754     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    755     REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence 
    756  
    757     REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s) 
    758     REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1) 
    759     REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm) 
    760     REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm) 
    761     REAL(rstd),INTENT(OUT) :: dps(iim*jjm) 
    762  
    763 ! temporary variable     
    764     INTEGER :: i,j,ij,l 
    765     REAL(rstd) :: p_ik, exner_ik 
    766     INTEGER    :: ij_omp_begin, ij_omp_end 
    767  
    768  
    769     CALL trace_start("compute_geopot") 
    770  
    771     CALL distrib_level(ij_end-ij_begin+1,ij_omp_begin,ij_omp_end) 
    772     ij_omp_begin=ij_omp_begin+ij_begin-1 
    773     ij_omp_end=ij_omp_end+ij_begin-1 
    774  
    775 !    REAL(rstd) :: wwuu(iim*3*jjm,llm+1) ! tmp var, don't know why but gain 30% on the whole code in opemp 
    776                                         ! need to be understood 
    777  
    778 !    wwuu=wwuu_out 
    779   CALL trace_start("compute_caldyn_vert") 
    780  
    781 !$OMP BARRIER    
    782 !!! cumulate mass flux convergence from top to bottom 
    783 !  IF (is_omp_level_master) THEN 
    784     DO  l = llm-1, 1, -1 
    785 !    IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    786  
    787 !!$OMP DO SCHEDULE(STATIC)  
    788 !DIR$ SIMD 
    789       DO ij=ij_omp_begin,ij_omp_end          
    790           convm(ij,l) = convm(ij,l) + convm(ij,l+1) 
    791       ENDDO 
    792     ENDDO 
    793 !  ENDIF 
    794  
    795 !$OMP BARRIER 
    796 ! FLUSH on convm 
    797 !!!!!!!!!!!!!!!!!!!!!!!!! 
    798  
    799 ! compute dps 
    800   IF (is_omp_first_level) THEN 
    801 !DIR$ SIMD 
    802      DO ij=ij_begin,ij_end          
    803         ! dps/dt = -int(div flux)dz 
    804         dps(ij) = convm(ij,1) * g  
    805     ENDDO 
    806   ENDIF 
    807    
    808 !!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC) 
    809   DO l=ll_beginp1,ll_end 
    810 !    IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    811 !DIR$ SIMD 
    812       DO ij=ij_begin,ij_end          
    813         ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 
    814         ! => w>0 for upward transport 
    815         wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l )  
    816     ENDDO 
    817   ENDDO 
    818  
    819  !--> flush wflux   
    820  !$OMP BARRIER 
    821  
    822   DO l=ll_begin,ll_endm1 
    823 !DIR$ SIMD 
    824      DO ij=ij_begin,ij_end          
    825         dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1))) 
    826     ENDDO 
    827   ENDDO 
    828    
    829   DO l=ll_beginp1,ll_end 
    830 !DIR$ SIMD 
    831       DO ij=ij_begin,ij_end          
    832         dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) ) 
    833     ENDDO 
    834   ENDDO 
    835  
    836    
    837 ! Compute vertical transport 
    838   DO l=ll_beginp1,ll_end 
    839 !DIR$ SIMD 
    840       DO ij=ij_begin,ij_end          
    841         wwuu(ij+u_right,l) = 0.5*( wflux(ij,l) + wflux(ij+t_right,l)) * (u(ij+u_right,l) - u(ij+u_right,l-1)) 
    842         wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1)) 
    843         wwuu(ij+u_ldown,l) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1)) 
    844      ENDDO 
    845    ENDDO 
    846  
    847  !--> flush wwuu   
    848  !$OMP BARRIER 
    849  
    850 ! Add vertical transport to du 
    851   DO l=ll_begin,ll_end 
    852 !DIR$ SIMD 
    853      DO ij=ij_begin,ij_end          
    854         du(ij+u_right, l )   = du(ij+u_right,l)  - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l)) 
    855         du(ij+u_lup, l )     = du(ij+u_lup,l)    - (wwuu(ij+u_lup,l+1)  + wwuu(ij+u_lup,l))   / (rhodz(ij,l)+rhodz(ij+t_lup,l)) 
    856         du(ij+u_ldown, l )   = du(ij+u_ldown,l)  - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l)) 
    857     ENDDO 
    858   ENDDO       
    859  
    860 !  DO l=ll_beginp1,ll_end 
    861 !!DIR$ SIMD 
    862 !      DO ij=ij_begin,ij_end          
    863 !        wwuu_out(ij+u_right,l) = wwuu(ij+u_right,l) 
    864 !        wwuu_out(ij+u_lup,l)   = wwuu(ij+u_lup,l) 
    865 !        wwuu_out(ij+u_ldown,l) = wwuu(ij+u_ldown,l) 
    866 !     ENDDO 
    867 !   ENDDO 
    868    
    869   CALL trace_end("compute_caldyn_vert") 
    870  
    871   END SUBROUTINE compute_caldyn_vert 
    872  
    873275!-------------------------------- Diagnostics ---------------------------- 
    874276 

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

    r360 r361  
    4848    USE time_mod 
    4949    USE disvert_mod 
    50     USE caldyn_mod 
     50    USE caldyn_hevi_mod 
    5151    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time 
    5252    INTEGER :: it,j,l, ind 
    … …  
    5454 
    5555    DO j=1,nb_stage 
    56        !       CALL caldyn_hevi((j==1) .AND. (MOD(it,itau_out)==0), & 
    57        !            taujj(j), f_phis, f_geopot, & 
    58        !            f_ps,f_mass,f_theta_rhodz,f_u, & 
    59        !            f_hflux, f_wflux, & 
    60        !            f_dps_slow(:,j), f_dmass_slow(:,j), f_dtheta_rhodz_slow(:,j), & 
    61        !            f_du_slow(:,j), f_du_fast(:,j) ) 
    62  
    63        CALL caldyn((j==1) .AND. (MOD(it,itau_out)==0), & 
     56       CALL caldyn_hevi((j==1) .AND. (MOD(it,itau_out)==0), taujj(j), & 
    6457            f_phis, f_ps,f_mass,f_theta_rhodz,f_u,f_q, & 
    6558            f_geopot, f_hflux, f_wflux, & 
    6659            f_dps_slow(:,j), f_dmass_slow(:,j), f_dtheta_rhodz_slow(:,j), & 
    67             f_du_slow(:,j)) 
    68  
    69        ! cumulate mass fluxes for transport scheme 
     60            f_du_slow(:,j), f_du_fast(:,j) ) 
     61       ! accumulate mass fluxes for transport scheme 
    7062       DO ind=1,ndomain 
    7163          IF (.NOT. assigned_domain(ind)) CYCLE 
    … …  
    7567          CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, wj(j), fluxt_zero(ind)) 
    7668       END DO 
    77  
    7869       ! update model state 
    7970       DO l=1,j 
    … …  
    8576          CALL update(bjl(l,j), f_theta_rhodz, f_dtheta_rhodz_slow(:,l)) 
    8677          CALL update(bjl(l,j), f_u, f_du_slow(:,l)) 
    87 !          CALL update(cjl(l,j), f_u, f_du_fast(:,l)) ! FIXME - only slow terms right now 
     78          CALL update(cjl(l,j), f_u, f_du_fast(:,l)) 
    8879       END DO 
    8980    END DO 
Note: See TracChangeset for help on using the changeset viewer.