Changeset 529 for codes/icosagcm/devel


Ignore:
Timestamp:
01/27/17 16:54:36 (7 years ago)
Author:
dubos
Message:

Updated devel to r526

Location:
codes/icosagcm/devel/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

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

    r413 r529  
    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/devel/src/caldyn_hevi.f90

    r404 r529  
    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/devel/src/caldyn_kernels.f90

    r428 r529  
    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             m = mass_dak(l)+(ps(ij)*g+ptop)*mass_dbk(l) ! ps is actually Ms 
     63             rhodz(ij,l) = m/g 
    6464             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    6565          ENDDO 
    … …  
    6767    ELSE ! Compute only theta 
    6868       DO l = ll_begin,ll_end 
    69           CALL test_message(req_u)  
     69 !         CALL test_message(req_u)  
    7070          !DIR$ SIMD 
    7171          DO ij=ij_begin_ext,ij_end_ext 
    … …  
    7575    END IF 
    7676 
    77     CALL wait_message(req_u) ! COM02 
     77 !   CALL wait_message(req_u) ! COM02 
    7878 
    7979!!! Compute shallow-water potential vorticity 
    … …  
    8181       !DIR$ SIMD 
    8282       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)     & 
     83          etav= 1./Av(ij+z_up)*(  ne_rup  * u(ij+u_rup,l)    & 
     84               + ne_left * u(ij+t_rup+u_left,l)              & 
     85               - ne_lup  * u(ij+u_lup,l) )                                
     86          hv =   Riv2(ij,vup)          * rhodz(ij,l)         & 
     87               + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)   & 
    8988               + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 
    90  
    9189          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)          & 
     90           
     91          etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)  & 
     92               + ne_right * u(ij+t_ldown+u_right,l)              & 
     93               - ne_rdown * u(ij+u_rdown,l) ) 
     94          hv =   Riv2(ij,vdown)        * rhodz(ij,l)          & 
    9895               + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  & 
    9996               + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 
    100  
    10197          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 
    102  
    10398       ENDDO 
    10499 
    … …  
    115110  END SUBROUTINE compute_pvort 
    116111 
    117   !************************* caldyn_horiz = caldyn_fast + caldyn_slow ********************** 
     112  !************** caldyn_horiz = caldyn_fast + caldyn_slow + caldyn_coriolis *************** 
    118113 
    119114  SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du) 
    … …  
    139134    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux 
    140135    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
     136    REAL(rstd) :: uu_right, uu_lup, uu_ldown 
    141137 
    142138    INTEGER :: i,j,ij,l 
    … …  
    149145    DO l = ll_begin, ll_end 
    150146!!!  Compute mass and theta fluxes 
    151        IF (caldyn_conserv==energy) CALL test_message(req_qu)  
     147!       IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    152148       !DIR$ SIMD 
    153149       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  
     150          uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) 
     151          uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) 
     152          uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) 
     153          uu_right= uu_right*le_de(ij+u_right) 
     154          uu_lup  = uu_lup  *le_de(ij+u_lup) 
     155          uu_ldown= uu_ldown*le_de(ij+u_ldown) 
     156          hflux(ij+u_right,l)=uu_right 
     157          hflux(ij+u_lup,l)  =uu_lup 
     158          hflux(ij+u_ldown,l)=uu_ldown 
     159!          hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 
     160!          hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 
     161!          hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 
    158162          Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l) 
    159163          Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l) 
    … …  
    189193    CASE(energy) ! energy-conserving TRiSK 
    190194 
    191        CALL wait_message(req_qu) ! COM03 
     195!       CALL wait_message(req_qu) ! COM03 
    192196 
    193197       DO l=ll_begin,ll_end 
    … …  
    205209                  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))+             & 
    206210                  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) 
     211             du(ij+u_right,l) = .5*uu 
    208212 
    209213             uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        & 
    … …  
    217221                  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)) +            & 
    218222                  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  
     223             du(ij+u_lup,l) = .5*uu 
    221224 
    222225             uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      & 
    … …  
    230233                  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)) +      & 
    231234                  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) 
     235             du(ij+u_ldown,l) = .5*uu 
    233236 
    234237          ENDDO 
    … …  
    251254                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   & 
    252255                  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) 
     256             du(ij+u_right,l) = qu(ij+u_right,l)*uu 
    254257 
    255258 
    … …  
    264267                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  & 
    265268                  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) 
     269             du(ij+u_lup,l) = qu(ij+u_lup,l)*uu 
    267270 
    268271             uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      & 
    … …  
    276279                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              & 
    277280                  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) 
     281             du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu 
    279282 
    280283          ENDDO 
    … …  
    290293          !DIR$ SIMD 
    291294          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 )   
     295             berni(ij,l) = pk(ij,l) + 1/(4*Ai(ij))*( & 
     296                  le_de(ij+u_right)*u(ij+u_right,l)**2 +    & 
     297                  le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
     298                  le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
     299                  le_de(ij+u_left)*u(ij+u_left,l)**2 +       & 
     300                  le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
     301                  le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    300302             ! from now on pk contains the vertically-averaged geopotential 
    301303             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 
    … …  
    308310          !DIR$ SIMD 
    309311          DO ij=ij_begin,ij_end          
    310  
    311312             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 )   
     313                  + 1/(4*Ai(ij))*( & 
     314                  le_de(ij+u_right)*u(ij+u_right,l)**2 +    & 
     315                  le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
     316                  le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
     317                  le_de(ij+u_left)*u(ij+u_left,l)**2 +       & 
     318                  le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
     319                  le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    318320          ENDDO 
    319321       ENDDO 
    … …  
    326328       DO ij=ij_begin,ij_end          
    327329 
    328           du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             & 
     330          du(ij+u_right,l) = du(ij+u_right,l) + (             & 
    329331               0.5*(theta(ij,l)+theta(ij+t_right,l))                & 
    330332               *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    & 
    … …  
    332334 
    333335 
    334           du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  & 
     336          du(ij+u_lup,l) = du(ij+u_lup,l) +  (                  & 
    335337               0.5*(theta(ij,l)+theta(ij+t_lup,l))                  & 
    336338               *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       & 
    337339               + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) ) 
    338340 
    339           du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             & 
     341          du(ij+u_ldown,l) = du(ij+u_ldown,l) + (             & 
    340342               0.5*(theta(ij,l)+theta(ij+t_ldown,l))                & 
    341343               *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     & 

  • codes/icosagcm/devel/src/caldyn_kernels_base.f90

    r479 r529  
    77  IMPLICIT NONE 
    88  PRIVATE 
    9  
    10   LOGICAL, PARAMETER, PUBLIC :: DEC = .TRUE. 
    119 
    1210  INTEGER, PARAMETER,PUBLIC :: energy=1, enstrophy=2 
    … …  
    219217       DO ij=ij_begin,ij_end          
    220218          ! dps/dt = -int(div flux)dz 
    221           IF(DEC) THEN 
    222              dps(ij) = convm(ij,1) 
    223           ELSE 
    224              dps(ij) = convm(ij,1) * g  
    225           END IF 
     219          dps(ij) = convm(ij,1) 
    226220       ENDDO 
    227221    ENDIF 

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

    r499 r529  
    3232          !DIR$ SIMD 
    3333          DO ij=ij_begin_ext,ij_end_ext 
    34              IF(DEC) THEN  ! ps is actually Ms 
    35                 m = mass_dak(l)+(ps(ij)*g+ptop)*mass_dbk(l) 
    36              ELSE 
    37                 m = mass_dak(l)+ps(ij)*mass_dbk(l) 
    38              END IF 
     34             m = mass_dak(l)+(ps(ij)*g+ptop)*mass_dbk(l) ! ps is actually Ms 
    3935             rhodz(ij,l) = m/g 
    4036          END DO 
    … …  
    6965       !DIR$ SIMD 
    7066       DO ij=ij_begin_ext,ij_end_ext 
    71           IF(DEC) THEN 
    72              etav= 1./Av(ij+z_up)*(  ne_rup  * u(ij+u_rup,l)         & 
    73                                    + ne_left * u(ij+t_rup+u_left,l)  & 
    74                                    - ne_lup  * u(ij+u_lup,l) )                                
    75  
    76              hv =   Riv2(ij,vup)          * rhodz(ij,l)           & 
    77                   + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     & 
    78                   + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 
    79              qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv 
    80  
    81              etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)          & 
    82                                       + ne_right * u(ij+t_ldown+u_right,l)  & 
    83                                       - ne_rdown * u(ij+u_rdown,l) ) 
    84              hv =   Riv2(ij,vdown)        * rhodz(ij,l)          & 
    85                   + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  & 
    86                   + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 
    87              qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 
    88           ELSE 
    89              etav= 1./Av(ij+z_up)*(  ne_rup  * u(ij+u_rup,l)        * de(ij+u_rup)         & 
    90                                    + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  & 
    91                                    - ne_lup  * u(ij+u_lup,l)        * de(ij+u_lup) )                                
    92  
    93              hv =   Riv2(ij,vup)          * rhodz(ij,l)           & 
    94                   + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     & 
    95                   + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 
    96              qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv 
    97  
    98              etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)          * de(ij+u_ldown)          & 
    99                                       + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  & 
    100                                       - ne_rdown * u(ij+u_rdown,l)          * de(ij+u_rdown) ) 
    101              hv =   Riv2(ij,vdown)        * rhodz(ij,l)          & 
    102                   + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  & 
    103                   + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 
    104              qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 
    105           END IF 
     67          etav= 1./Av(ij+z_up)*(  ne_rup  * u(ij+u_rup,l)   & 
     68                  + ne_left * u(ij+t_rup+u_left,l)          & 
     69                  - ne_lup  * u(ij+u_lup,l) )                                
     70          hv =   Riv2(ij,vup)          * rhodz(ij,l)        & 
     71               + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)  & 
     72               + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 
     73          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv 
     74           
     75          etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)   & 
     76               + ne_right * u(ij+t_ldown+u_right,l)               & 
     77               - ne_rdown * u(ij+u_rdown,l) ) 
     78          hv =   Riv2(ij,vdown)        * rhodz(ij,l)              & 
     79               + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)      & 
     80               + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 
     81          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 
    10682       ENDDO 
    10783 
    … …  
    649625  END SUBROUTINE compute_caldyn_Coriolis 
    650626 
    651   SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du) 
     627  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du, zero) 
     628    LOGICAL, INTENT(IN) :: zero 
    652629    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity" 
    653630    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    654631    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s 
    655     REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm) 
     632    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm) 
    656633     
    657634    REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function 
    … …  
    690667       ENDDO 
    691668       ! Compute du=-grad(Bernoulli) 
    692        !DIR$ SIMD 
    693        DO ij=ij_begin,ij_end 
     669       IF(zero) THEN 
     670          !DIR$ SIMD 
     671          DO ij=ij_begin,ij_end 
    694672             du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right)) 
    695673             du(ij+u_lup,l)   = ne_lup*(berni(ij)-berni(ij+t_lup)) 
    696674             du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown)) 
    697        END DO 
     675          END DO 
     676       ELSE 
     677          !DIR$ SIMD 
     678          DO ij=ij_begin,ij_end 
     679             du(ij+u_right,l) = du(ij+u_right,l) + & 
     680                  ne_right*(berni(ij)-berni(ij+t_right)) 
     681             du(ij+u_lup,l)   = du(ij+u_lup,l) + & 
     682                  ne_lup*(berni(ij)-berni(ij+t_lup)) 
     683             du(ij+u_ldown,l) = du(ij+u_ldown,l) + & 
     684                  ne_ldown*(berni(ij)-berni(ij+t_ldown)) 
     685          END DO 
     686       END IF 
    698687    END DO 
    699688 

  • codes/icosagcm/devel/src/dissip_gcm.f90

    r487 r529  
    77  TYPE(t_field),POINTER,SAVE :: f_due_diss2(:) 
    88 
    9   TYPE(t_field),POINTER,SAVE :: f_theta(:), f_phi(:), f_pk(:), f_pks(:), f_p(:) 
    109  TYPE(t_field),POINTER,SAVE :: f_dtheta_diss(:) 
    1110  TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz_diss(:) 
    … …  
    3534  INTEGER, SAVE :: rayleigh_friction_type, rayleigh_shear 
    3635!$OMP THREADPRIVATE(rayleigh_friction_type) 
     36!$OMP THREADPRIVATE(rayleigh_shear) 
    3737  REAL, SAVE    :: rayleigh_tau 
    38 !$OMP THREADPRIVATE(rayleigh_shear) 
     38!$OMP THREADPRIVATE(rayleigh_tau) 
    3939   
    4040  REAL,SAVE    :: dtdissip 
    … …  
    4949    CALL allocate_field(f_due_diss1,field_u,type_real,llm) 
    5050    CALL allocate_field(f_due_diss2,field_u,type_real,llm) 
    51     CALL allocate_field(f_theta,field_t,type_real,llm) 
    5251    CALL allocate_field(f_dtheta_diss,field_t,type_real,llm) 
    5352    CALL allocate_field(f_dtheta_rhodz_diss,field_t,type_real,llm) 
    54  
    55     CALL allocate_field(f_phi,field_t,type_real,llm) 
    56     CALL allocate_field(f_pk,field_t,type_real,llm) 
    57     CALL allocate_field(f_p,field_t,type_real,llm+1) 
    58     CALL allocate_field(f_pks,field_t,type_real) 
    59      
    6053    ALLOCATE(tau_graddiv(llm)) 
    6154    ALLOCATE(tau_gradrot(llm))     
    … …  
    435428      mintau=MIN(mintau,tau_divgrad(l)) 
    436429    ENDDO 
    437         
     430 
     431    IF(rayleigh_friction_type>0) mintau=MIN(mintau,rayleigh_tau) 
     432 
    438433    IF(mintau>0) THEN 
    439434       itau_dissip=INT(mintau/dt) 
    … …  
    444439    END IF 
    445440    itau_dissip=MAX(1,itau_dissip) 
    446     IF (is_master) PRINT *,"mintau ",mintau,"itau_dissip",itau_dissip," dtdissip ",dtdissip 
     441    IF (is_master) PRINT *,"rayleigh_tau",rayleigh_tau, "mintau ",mintau, & 
     442         "itau_dissip",itau_dissip," dtdissip ",dtdissip 
    447443 
    448444  END SUBROUTINE init_dissip  
    449445  
    450446   
    451   SUBROUTINE dissip(f_ue,f_due,f_mass,f_phis,f_theta_rhodz,f_dtheta_rhodz) 
     447  SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due) 
    452448  USE icosa 
    453449  USE theta2theta_rhodz_mod 
    … …  
    459455  USE omp_para 
    460456  IMPLICIT NONE 
    461     TYPE(t_field),POINTER :: f_ue(:) 
    462     TYPE(t_field),POINTER :: f_due(:) 
    463     TYPE(t_field),POINTER :: f_mass(:), f_phis(:) 
    464     TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    465     TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 
     457    TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:) 
     458    TYPE(t_field),POINTER :: f_theta_rhodz(:), f_dtheta_rhodz(:) 
     459    TYPE(t_field),POINTER :: f_ue(:), f_due(:) 
    466460 
    467461    REAL(rstd),POINTER         :: due(:,:) 
    … …  
    483477    CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss) 
    484478     
    485 ! later for openmp     
    486 !    IF(rayleigh_friction_type>0) THEN 
    487 !       CALL pression(f_ps, f_p) 
    488 !       CALL exner(f_ps, f_p, f_pks, f_pk) 
    489 !       CALL geopotential(f_phis,f_pks,f_pk,f_theta,f_phi) 
    490 !    END IF 
    491  
    492479    DO ind=1,ndomain 
    493480      IF (.NOT. assigned_domain(ind)) CYCLE 
    … …  
    515502!      due=0 
    516503 
    517 ! later for openmp   
    518 !      IF(rayleigh_friction_type>0) THEN 
    519 !         phi=f_phi(ind) 
    520 !         ue=f_ue(ind) 
    521 !         DO l=1,llm 
    522 !            DO j=jj_begin,jj_end 
    523 !               DO i=ii_begin,ii_end 
    524 !                  n=(j-1)*iim+i 
    525 !                  CALL relax(t_right, u_right) 
    526 !                  CALL relax(t_lup,   u_lup) 
    527 !                  CALL relax(t_ldown, u_ldown) 
    528 !               ENDDO 
    529 !            ENDDO 
    530 !         END DO 
    531 !      END IF 
     504      IF(rayleigh_friction_type>0) THEN 
     505         phi=f_geopot(ind) 
     506         ue=f_ue(ind) 
     507         DO l=ll_begin,ll_end 
     508            DO ij=ij_begin,ij_end 
     509               CALL relax(t_right, u_right) 
     510               CALL relax(t_lup,   u_lup) 
     511               CALL relax(t_ldown, u_ldown) 
     512            ENDDO 
     513         END DO 
     514      END IF 
    532515   END DO 
    533516 
    … …  
    558541           fz = sin((pi/2)*(z-zh)/(ztop-zh)) 
    559542           fz = fz*fz/rayleigh_tau 
    560 !           fz = 1/rayleigh_tau 
    561            due(nn,l) = due(nn,l) - fz*(ue(nn,l)-uref) 
     543           due(nn,l) = due(nn,l) - itau_dissip*fz*(ue(nn,l)-uref) 
     544!           fz = 1./rayleigh_tau 
    562545!           due(nn,l) = due(nn,l) - fz*ue(nn,l) 
    563546         END IF 

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

    r487 r529  
    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/devel/src/hevi_scheme.f90

    r480 r529  
    44  USE field_mod 
    55  USE euler_scheme_mod 
    6   USE caldyn_kernels_base_mod, ONLY : DEC 
    76  IMPLICIT NONE 
    87  PRIVATE 
    … …  
    6059    REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:) 
    6160 
    62     IF(DEC) CALL legacy_to_DEC(f_ps, f_u) 
     61    CALL legacy_to_DEC(f_ps, f_u) 
    6362    DO j=1,nb_stage 
    6463       CALL caldyn_hevi((j==1) .AND. (MOD(it,itau_out)==0), taujj(j), & 
    … …  
    9796       END DO 
    9897    END DO 
    99     IF(DEC) CALL DEC_to_legacy(f_ps, f_u) 
     98    CALL DEC_to_legacy(f_ps, f_u) 
    10099  END SUBROUTINE HEVI_scheme 
    101100   

  • codes/icosagcm/devel/src/observable.f90

    r482 r529  
    7777    END IF 
    7878 
    79     CALL divide_by_mass(1, f_mass, f_theta_rhodz, f_buf_i) 
    80     IF(init) THEN 
    81        CALL output_field("theta_init",f_buf_i) 
    82     ELSE 
    83        CALL output_field("theta",f_buf_i) 
    84     END IF 
    85  
    8679    IF(nqdyn>1) THEN 
    8780       CALL divide_by_mass(2, f_mass, f_theta_rhodz, f_buf_i) 
    … …  
    9386    END IF 
    9487 
    95 !    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) 
    96 !    CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
    97     CALL diagnose_temperature(f_ps, f_theta_rhodz, f_q, f_buf_i) 
     88    CALL divide_by_mass(1, f_mass, f_theta_rhodz, f_buf_i) 
     89    IF(init) THEN 
     90       CALL output_field("theta_init",f_buf_i) 
     91    ELSE 
     92       CALL output_field("theta",f_buf_i) 
     93    END IF 
     94 
     95    CALL pression_mid(f_ps, f_pmid) 
     96    CALL diagnose_temperature(f_pmid, f_q, f_buf_i) ! f_buf_i : IN = theta, out = T 
    9897 
    9998    IF(init) THEN 
    … …  
    112111    CALL transfert_request(f_buf_uh,req_e1_vect)  
    113112    CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat) 
    114     CALL pression_mid(f_ps, f_pmid) 
    115113    IF(init) THEN 
    116114       CALL output_field("uz_init",f_buf_i) 
    … …  
    236234  END SUBROUTINE compute_prognostic_vel_to_horiz 
    237235 
    238   SUBROUTINE diagnose_temperature(f_ps,f_theta_rhodz,f_q,f_temp) 
     236  SUBROUTINE diagnose_temperature(f_pmid,f_q,f_temp) 
    239237    USE icosa 
    240238    USE pression_mod 
    241239    IMPLICIT NONE 
    242     TYPE(t_field), POINTER :: f_ps(:)           ! IN 
    243     TYPE(t_field), POINTER :: f_theta_rhodz(:)  ! IN 
     240    TYPE(t_field), POINTER :: f_pmid(:)           ! IN 
    244241    TYPE(t_field), POINTER :: f_q(:)            ! IN 
    245     TYPE(t_field), POINTER :: f_temp(:)         ! OUT 
    246      
    247     REAL(rstd), POINTER :: ps(:) 
    248     REAL(rstd), POINTER :: theta_rhodz(:,:,:) 
     242    TYPE(t_field), POINTER :: f_temp(:)         ! INOUT 
     243     
     244    REAL(rstd), POINTER :: pmid(:,:) 
    249245    REAL(rstd), POINTER :: q(:,:,:) 
    250246    REAL(rstd), POINTER :: temp(:,:) 
    … …  
    255251       CALL swap_dimensions(ind) 
    256252       CALL swap_geometry(ind) 
    257        ps=f_ps(ind) 
    258        theta_rhodz=f_theta_rhodz(ind) 
     253       pmid=f_pmid(ind) 
    259254       q=f_q(ind) 
    260255       temp=f_temp(ind) 
    261        CALL compute_diagnose_temp(ps,theta_rhodz,q,temp) 
     256       CALL compute_diagnose_temp(pmid,q,temp) 
    262257    END DO 
    263258  END SUBROUTINE diagnose_temperature 
    264259   
    265   SUBROUTINE compute_diagnose_temp(ps,theta_rhodz,q,temp) 
     260  SUBROUTINE compute_diagnose_temp(pmid,q,temp) 
    266261    USE omp_para 
    267262    USE pression_mod 
    268     REAL(rstd),INTENT(IN)  :: ps(iim*jjm) 
    269     REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm,nqdyn) 
    270     REAL(rstd),INTENT(IN)  :: q(iim*jjm,llm,nqtot) 
    271     REAL(rstd),INTENT(OUT) :: temp(iim*jjm,llm) 
    272  
    273     REAL(rstd)  :: p(iim*jjm,llm+1) 
     263    REAL(rstd),INTENT(IN)    :: pmid(iim*jjm,llm) 
     264    REAL(rstd),INTENT(IN)    :: q(iim*jjm,llm,nqtot) 
     265    REAL(rstd),INTENT(INOUT) :: temp(iim*jjm,llm) 
     266 
    274267    REAL(rstd) :: Rd, p_ik, theta_ik, temp_ik, qv, chi, Rmix 
    275268    INTEGER :: ij,l 
    276269 
    277270    Rd = kappa*cpp 
    278     CALL compute_pression(ps,p,0) 
    279271    DO l=ll_begin,ll_end 
    280272       DO ij=ij_begin,ij_end 
    281           p_ik = .5*(p(ij,l)+p(ij,l+1)) 
    282           theta_ik = g*theta_rhodz(ij,l,1)/(p(ij,l)-p(ij,l+1)) 
     273          p_ik = pmid(ij,l) 
     274          theta_ik = temp(ij,l) 
    283275          qv = q(ij,l,1) ! water vaper mixing ratio = mv/md 
    284276          SELECT CASE(caldyn_thermo) 

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

    r528 r529  
    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 
    … …  
    284284               f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 
    285285        
    286           CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz) 
     286          CALL dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_u, f_dtheta_rhodz,f_du) 
    287287           
    288288          CALL euler_scheme(.FALSE.)  ! update only u, theta 
Note: See TracChangeset for help on using the changeset viewer.