Changeset 519 for codes


Ignore:
Timestamp:
12/24/16 02:33:07 (7 years ago)
Author:
dubos
Message:

Fixed RK2.5 - cleanup to follow

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

Legend:

Unmodified
Added
Removed

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

    r413 r519  
    22  USE icosa 
    33  USE transfert_mod 
     4  USE caldyn_kernels_hevi_mod 
    45  USE caldyn_kernels_base_mod 
    56  USE caldyn_kernels_mod 
    … …  
    190191 
    191192! temporary shared variable 
    192     REAL(rstd),POINTER  :: theta(:,:)   
     193    REAL(rstd),POINTER  :: theta(:,:,:)   
    193194    REAL(rstd),POINTER  :: pk(:,:) 
    194195    REAL(rstd),POINTER  :: geopot(:,:) 
    … …  
    225226    IF(caldyn_eta==eta_mass) THEN 
    226227       CALL send_message(f_ps,req_ps) ! COM00 
     228       CALL wait_message(req_ps) ! COM00 
    227229    ELSE 
    228230       CALL send_message(f_mass,req_mass) ! COM00 
     231       CALL wait_message(req_mass) ! COM00 
    229232    END IF 
    230233     
    231234    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01 
     235    CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort 
    232236    CALL send_message(f_u,req_u) ! COM02 
     237    CALL wait_message(req_u) ! COM02 
     238 
     239    IF(.NOT.hydrostatic) THEN 
     240       STOP 'caldyn_gcm may not be used yet when non-hydrostatic' 
     241    END IF 
    233242 
    234243    SELECT CASE(caldyn_conserv) 
    … …  
    245254          qu=f_qu(ind) 
    246255          qv=f_qv(ind) 
     256          pk = f_pk(ind) 
     257          geopot = f_geopot(ind)   
     258          hflux=f_hflux(ind) 
     259          convm = f_dmass(ind) 
     260          dtheta_rhodz=f_dtheta_rhodz(ind) 
     261          du=f_du(ind) 
    247262          CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv) ! COM00 COM01 COM02 
    248        ENDDO 
    249 !       CALL checksum(f_mass) 
    250 !       CALL checksum(f_theta) 
    251  
    252        CALL send_message(f_qu,req_qu) ! COM03 wait_message in caldyn_horiz 
    253 !       CALL wait_message(req_qu) 
     263!          CALL compute_theta(ps,theta_rhodz, mass,theta) 
     264!          CALL compute_pvort_only(u,mass,qu,qv) 
     265 
     266          CALL compute_geopot(mass,theta, ps,pk,geopot) 
     267!          du(:,:)=0. 
     268!          CALL compute_caldyn_fast(0.,u,mass,theta,pk,geopot,du) 
     269       ENDDO        
     270 
     271       CALL send_message(f_u,req_u) ! COM02 
     272       CALL wait_message(req_u) ! COM02 
     273       CALL send_message(f_qu,req_qu) ! COM03 
     274       CALL wait_message(req_qu) ! COM03 
    254275 
    255276       DO ind=1,ndomain 
    … …  
    259280          ps=f_ps(ind) 
    260281          u=f_u(ind) 
     282          theta_rhodz = f_theta_rhodz(ind) 
    261283          mass=f_mass(ind) 
    262284          theta = f_theta(ind) 
    263285          qu=f_qu(ind) 
     286          qv=f_qv(ind) 
    264287          pk = f_pk(ind) 
    265288          geopot = f_geopot(ind)   
    266           CALL compute_geopot(ps,mass,theta, pk,geopot) 
    267289          hflux=f_hflux(ind) 
    268290          convm = f_dmass(ind) 
    269291          dtheta_rhodz=f_dtheta_rhodz(ind) 
    270292          du=f_du(ind) 
     293 
    271294          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du) 
     295!          CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .FALSE.) ! FIXME 
     296!          CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)  
    272297          IF(caldyn_eta==eta_mass) THEN 
    273298             wflux=f_wflux(ind) 
    … …  
    278303       ENDDO        
    279304    
    280 !       CALL checksum(f_geopot) 
    281 !       CALL checksum(f_dmass) 
    282 !       CALL checksum(f_pk) 
    283 !       CALL checksum(f_pk) 
    284         
    285305    CASE(enstrophy) ! enstrophy-conserving 
    286306       DO ind=1,ndomain 

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

    r404 r519  
    155155 
    156156       IF(hydrostatic) THEN 
    157           CALL compute_caldyn_slow_hydro(u,mass,hflux,du) 
     157          CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .TRUE.) 
    158158       ELSE 
    159159          W = f_W(ind) 

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

    r428 r519  
    3232  SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv) 
    3333    USE icosa 
    34     USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass 
     34    USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop 
    3535    USE trace 
    3636    USE omp_para 
    … …  
    4848    CALL trace_start("compute_pvort")   
    4949 
    50     IF(caldyn_eta==eta_mass) THEN 
    51        CALL wait_message(req_ps)  ! COM00 
    52     ELSE 
    53        CALL wait_message(req_mass) ! COM00 
    54     END IF 
    55     CALL wait_message(req_theta_rhodz) ! COM01 
     50!    IF(caldyn_eta==eta_mass) THEN 
     51!       CALL wait_message(req_ps)  ! COM00 
     52!    ELSE 
     53!       CALL wait_message(req_mass) ! COM00 
     54!    END IF 
     55!    CALL wait_message(req_theta_rhodz) ! COM01 
    5656 
    5757    IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta 
    5858       DO l = ll_begin,ll_end 
    59           CALL test_message(req_u)  
     59!          CALL test_message(req_u)  
    6060          !DIR$ SIMD 
    6161          DO ij=ij_begin_ext,ij_end_ext 
    62              m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 
    63              rhodz(ij,l) = m 
     62             IF(DEC) THEN  ! ps is actually Ms 
     63                m = mass_dak(l)+(ps(ij)*g+ptop)*mass_dbk(l) 
     64             ELSE 
     65                m = mass_dak(l)+ps(ij)*mass_dbk(l) 
     66             END IF 
     67             rhodz(ij,l) = m/g 
    6468             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    6569          ENDDO 
    … …  
    6771    ELSE ! Compute only theta 
    6872       DO l = ll_begin,ll_end 
    69           CALL test_message(req_u)  
     73 !         CALL test_message(req_u)  
    7074          !DIR$ SIMD 
    7175          DO ij=ij_begin_ext,ij_end_ext 
    … …  
    7579    END IF 
    7680 
    77     CALL wait_message(req_u) ! COM02 
     81 !   CALL wait_message(req_u) ! COM02 
    7882 
    7983!!! Compute shallow-water potential vorticity 
    … …  
    8185       !DIR$ SIMD 
    8286       DO ij=ij_begin_ext,ij_end_ext 
    83           etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         & 
    84                + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  & 
    85                - ne_lup        * u(ij+u_lup,l)        * de(ij+u_lup) )                                
    86  
    87           hv =  Riv2(ij,vup)          * rhodz(ij,l)            & 
    88                + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     & 
    89                + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 
    90  
    91           qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv 
    92  
    93           etav = 1./Av(ij+z_down)*(  ne_ldown         * u(ij+u_ldown,l)          * de(ij+u_ldown)          & 
    94                + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  & 
    95                - ne_rdown         * u(ij+u_rdown,l)          * de(ij+u_rdown) ) 
    96  
    97           hv =  Riv2(ij,vdown)        * rhodz(ij,l)          & 
    98                + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  & 
    99                + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 
    100  
    101           qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 
    102  
     87          IF(DEC) THEN 
     88             etav= 1./Av(ij+z_up)*(  ne_rup  * u(ij+u_rup,l)         & 
     89                                   + ne_left * u(ij+t_rup+u_left,l)  & 
     90                                   - ne_lup  * u(ij+u_lup,l) )                                
     91 
     92             hv =   Riv2(ij,vup)          * rhodz(ij,l)           & 
     93                  + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     & 
     94                  + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 
     95             qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv 
     96 
     97             etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)          & 
     98                                      + ne_right * u(ij+t_ldown+u_right,l)  & 
     99                                      - ne_rdown * u(ij+u_rdown,l) ) 
     100             hv =   Riv2(ij,vdown)        * rhodz(ij,l)          & 
     101                  + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  & 
     102                  + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 
     103             qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 
     104          ELSE 
     105             etav= 1./Av(ij+z_up)*(  ne_rup  * u(ij+u_rup,l)        * de(ij+u_rup)         & 
     106                                   + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  & 
     107                                   - ne_lup  * u(ij+u_lup,l)        * de(ij+u_lup) )                                
     108 
     109             hv =   Riv2(ij,vup)          * rhodz(ij,l)           & 
     110                  + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     & 
     111                  + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 
     112             qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv 
     113 
     114             etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)          * de(ij+u_ldown)          & 
     115                                      + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  & 
     116                                      - ne_rdown * u(ij+u_rdown,l)          * de(ij+u_rdown) ) 
     117             hv =   Riv2(ij,vdown)        * rhodz(ij,l)          & 
     118                  + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  & 
     119                  + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 
     120             qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 
     121          END IF 
    103122       ENDDO 
    104123 
    … …  
    115134  END SUBROUTINE compute_pvort 
    116135 
    117   !************************* caldyn_horiz = caldyn_fast + caldyn_slow ********************** 
     136  !************** caldyn_horiz = caldyn_fast + caldyn_slow + caldyn_coriolis *************** 
    118137 
    119138  SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du) 
    … …  
    139158    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux 
    140159    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
     160    REAL(rstd) :: uu_right, uu_lup, uu_ldown 
    141161 
    142162    INTEGER :: i,j,ij,l 
    … …  
    149169    DO l = ll_begin, ll_end 
    150170!!!  Compute mass and theta fluxes 
    151        IF (caldyn_conserv==energy) CALL test_message(req_qu)  
     171!       IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    152172       !DIR$ SIMD 
    153173       DO ij=ij_begin_ext,ij_end_ext          
    154           hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 
    155           hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 
    156           hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 
    157  
     174          uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) 
     175          uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) 
     176          uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) 
     177          uu_right= uu_right*le_de(ij+u_right) 
     178          uu_lup  = uu_lup  *le_de(ij+u_lup) 
     179          uu_ldown= uu_ldown*le_de(ij+u_ldown) 
     180          hflux(ij+u_right,l)=uu_right 
     181          hflux(ij+u_lup,l)  =uu_lup 
     182          hflux(ij+u_ldown,l)=uu_ldown 
     183!          hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 
     184!          hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 
     185!          hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 
    158186          Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l) 
    159187          Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l) 
    … …  
    189217    CASE(energy) ! energy-conserving TRiSK 
    190218 
    191        CALL wait_message(req_qu) ! COM03 
     219!       CALL wait_message(req_qu) ! COM03 
    192220 
    193221       DO l=ll_begin,ll_end 
    … …  
    205233                  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))+             & 
    206234                  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)) 
    207              du(ij+u_right,l) = .5*uu/de(ij+u_right) 
     235             du(ij+u_right,l) = .5*uu 
    208236 
    209237             uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        & 
    … …  
    217245                  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)) +            & 
    218246                  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)) 
    219              du(ij+u_lup,l) = .5*uu/de(ij+u_lup) 
    220  
     247             du(ij+u_lup,l) = .5*uu 
    221248 
    222249             uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      & 
    … …  
    230257                  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)) +      & 
    231258                  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)) 
    232              du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown) 
     259             du(ij+u_ldown,l) = .5*uu 
    233260 
    234261          ENDDO 
    … …  
    251278                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   & 
    252279                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l) 
    253              du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right) 
     280             du(ij+u_right,l) = qu(ij+u_right,l)*uu 
    254281 
    255282 
    … …  
    264291                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  & 
    265292                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l) 
    266              du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup) 
     293             du(ij+u_lup,l) = qu(ij+u_lup,l)*uu 
    267294 
    268295             uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      & 
    … …  
    276303                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              & 
    277304                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l) 
    278              du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown) 
     305             du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu 
    279306 
    280307          ENDDO 
    … …  
    290317          !DIR$ SIMD 
    291318          DO ij=ij_begin,ij_end          
    292  
    293              berni(ij,l) = pk(ij,l) + & 
    294                   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    & 
    295                   le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
    296                   le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
    297                   le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       & 
    298                   le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
    299                   le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
     319             berni(ij,l) = pk(ij,l) + 1/(4*Ai(ij))*( & 
     320                  le_de(ij+u_right)*u(ij+u_right,l)**2 +    & 
     321                  le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
     322                  le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
     323                  le_de(ij+u_left)*u(ij+u_left,l)**2 +       & 
     324                  le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
     325                  le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    300326             ! from now on pk contains the vertically-averaged geopotential 
    301327             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 
    … …  
    308334          !DIR$ SIMD 
    309335          DO ij=ij_begin,ij_end          
    310  
    311336             berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
    312                   + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    & 
    313                   le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
    314                   le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
    315                   le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       & 
    316                   le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
    317                   le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
     337                  + 1/(4*Ai(ij))*( & 
     338                  le_de(ij+u_right)*u(ij+u_right,l)**2 +    & 
     339                  le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
     340                  le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
     341                  le_de(ij+u_left)*u(ij+u_left,l)**2 +       & 
     342                  le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
     343                  le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    318344          ENDDO 
    319345       ENDDO 
    … …  
    326352       DO ij=ij_begin,ij_end          
    327353 
    328           du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             & 
     354          du(ij+u_right,l) = du(ij+u_right,l) + (             & 
    329355               0.5*(theta(ij,l)+theta(ij+t_right,l))                & 
    330356               *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    & 
    … …  
    332358 
    333359 
    334           du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  & 
     360          du(ij+u_lup,l) = du(ij+u_lup,l) +  (                  & 
    335361               0.5*(theta(ij,l)+theta(ij+t_lup,l))                  & 
    336362               *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       & 
    337363               + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) ) 
    338364 
    339           du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             & 
     365          du(ij+u_ldown,l) = du(ij+u_ldown,l) + (             & 
    340366               0.5*(theta(ij,l)+theta(ij+t_ldown,l))                & 
    341367               *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     & 

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

    r499 r519  
    649649  END SUBROUTINE compute_caldyn_Coriolis 
    650650 
    651   SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du) 
     651  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du, zero) 
     652    LOGICAL, INTENT(IN) :: zero 
    652653    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity" 
    653654    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    654655    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s 
    655     REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm) 
     656    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm) 
    656657     
    657658    REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function 
    … …  
    690691       ENDDO 
    691692       ! Compute du=-grad(Bernoulli) 
    692        !DIR$ SIMD 
    693        DO ij=ij_begin,ij_end 
     693       IF(zero) THEN 
     694          !DIR$ SIMD 
     695          DO ij=ij_begin,ij_end 
    694696             du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right)) 
    695697             du(ij+u_lup,l)   = ne_lup*(berni(ij)-berni(ij+t_lup)) 
    696698             du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown)) 
    697        END DO 
     699          END DO 
     700       ELSE 
     701          !DIR$ SIMD 
     702          DO ij=ij_begin,ij_end 
     703             du(ij+u_right,l) = du(ij+u_right,l) + & 
     704                  ne_right*(berni(ij)-berni(ij+t_right)) 
     705             du(ij+u_lup,l)   = du(ij+u_lup,l) + & 
     706                  ne_lup*(berni(ij)-berni(ij+t_lup)) 
     707             du(ij+u_ldown,l) = du(ij+u_ldown,l) + & 
     708                  ne_ldown*(berni(ij)-berni(ij+t_ldown)) 
     709          END DO 
     710       END IF 
    698711    END DO 
    699712 

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

    r487 r519  
    2626    REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:), du(:,:) 
    2727    REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), massm1(:,:), massm2(:,:), dmass(:,:) 
    28     REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:), dtheta_rhodz(:,:) 
     28    REAL(rstd),POINTER :: theta_rhodz(:,:,:) , theta_rhodzm1(:,:,:), theta_rhodzm2(:,:,:), dtheta_rhodz(:,:,:) 
    2929    REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:) 
    3030    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time 
    3131    INTEGER :: it,stage 
     32     
     33    CALL legacy_to_DEC(f_ps, f_u) 
    3234    DO stage=1,nb_stage 
    3335       !       CALL checksum(f_ps) 
    … …  
    5355       END SELECT 
    5456    END DO 
     57    CALL DEC_to_legacy(f_ps, f_u) 
    5558 
    5659  CONTAINS 
    … …  
    142145                  um1(ij+u_lup,l)=u(ij+u_lup,l) 
    143146                  um1(ij+u_ldown,l)=u(ij+u_ldown,l) 
    144                   theta_rhodzm1(ij,l)=theta_rhodz(ij,l) 
     147                  theta_rhodzm1(ij,l,1)=theta_rhodz(ij,l,1) 
    145148               ENDDO 
    146149            ENDDO 
    … …  
    153156               u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l) 
    154157               u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l) 
    155                theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l) 
     158               theta_rhodz(ij,l,1)=theta_rhodzm1(ij,l,1)+tau*dtheta_rhodz(ij,l,1) 
    156159            ENDDO 
    157160         ENDDO 
    … …  
    188191 
    189192         IF (stage==1) THEN 
    190             psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) 
     193            psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:,:)=theta_rhodz(:,:,:) 
    191194            ps(:)=psm1(:)+tau*dps(:) 
    192195            u(:,:)=um1(:,:)+tau*du(:,:) 
    193             theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) 
     196            theta_rhodz(:,:,:)=theta_rhodzm1(:,:,:)+tau*dtheta_rhodz(:,:,:) 
    194197 
    195198         ELSE IF (stage==2) THEN 
    … …  
    197200            ps(:)=psm1(:)+tau*dps(:) 
    198201            u(:,:)=um1(:,:)+tau*du(:,:) 
    199             theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) 
    200  
    201             psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)  
    202             psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:)  
     202            theta_rhodz(:,:,:)=theta_rhodzm1(:,:,:)+tau*dtheta_rhodz(:,:,:) 
     203 
     204            psm2(:)=psm1(:) ; theta_rhodzm2(:,:,:)=theta_rhodzm1(:,:,:) ; um2(:,:)=um1(:,:)  
     205            psm1(:)=ps(:) ; theta_rhodzm1(:,:,:)=theta_rhodz(:,:,:) ; um1(:,:)=u(:,:)  
    203206 
    204207         ELSE  
    … …  
    206209            ps(:)=psm2(:)+2*tau*dps(:) 
    207210            u(:,:)=um2(:,:)+2*tau*du(:,:) 
    208             theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:) 
    209  
    210             psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)  
    211             psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:)  
     211            theta_rhodz(:,:,:)=theta_rhodzm2(:,:,:)+2*tau*dtheta_rhodz(:,:,:) 
     212 
     213            psm2(:)=psm1(:) ; theta_rhodzm2(:,:,:)=theta_rhodzm1(:,:,:) ; um2(:,:)=um1(:,:)  
     214            psm1(:)=ps(:) ; theta_rhodzm1(:,:,:)=theta_rhodz(:,:,:) ; um1(:,:)=u(:,:)  
    212215 
    213216         ENDIF 

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

    r516 r519  
    211211    IF(positive_theta) CALL copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q) 
    212212 
     213    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0)   
     214 
     215    CALL trace_on 
     216 
     217    IF (xios_output) THEN ! we must call update_calendar before any XIOS output 
     218      IF (is_omp_master) CALL xios_update_calendar(1) 
     219    END IF 
     220     
     221    IF (write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q)  ! FIXME : write_start undefined 
     222    
     223    CALL write_output_fields_basic(.TRUE., f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q) 
     224 
    213225    !$OMP MASTER 
    214226    CALL SYSTEM_CLOCK(start_clock) 
    215227    CALL SYSTEM_CLOCK(count_rate=rate_clock) 
    216228    !$OMP END MASTER    
    217  
    218     CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0)   
    219  
    220     CALL trace_on 
    221  
    222     IF (xios_output) THEN ! we must call update_calendar before any XIOS output 
    223       IF (is_omp_master) CALL xios_update_calendar(1) 
    224     END IF 
    225      
    226     IF (write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q)  ! FIXME : write_start undefined 
    227     
    228     CALL write_output_fields_basic(.TRUE., f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q) 
    229229 
    230230    DO it=itau0+1,itau0+itaumax 
Note: See TracChangeset for help on using the changeset viewer.